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ABSTRACT 


The objective of this research is to advance the development of a micro-analysis 
technique for characterizing the charge transport properties in bulk semiconductor 
materials for room temperature nuclear radiation detection. The technique is applied to 
bulk semi-insulating thallium bromide (TIBr) and cadmium zinc telluride (CZT) and uses 
a two-dimensional diffusion model and transport imaging to make rapid contact-free 
measurements of the magnitude and spatial variation in the mobility-lifetime {pi z) product 
at 2 pm resolution. The pi z product is a key measure of charge transport, and a uniform 
pi z product is critical for optimum energy resolution in gamma ray detectors. Spatial 
variations in the ambipolar diffusion length are observed in TIBr on a scale of ~10 pm, 
and the cathodoluminescence (CL) spectrum in TIBr at 5 K is reported. Using CL and 
photo-induced conductivity transient spectroscopy (PICTS), an empirical energy level 
diagram for Se and Pb defect levels in TIBr is reported. The pi z product in TIBr is 
generally found to decrease with increasing impurities/defects. The pi z product in TIBr 
and CZT is found to decrease with increasing temperature over the ranges of 8 K-102 K 
and 5 K-60 K, respectively. Transport imaging is applied for the first time to investigate 
the effects of Te inclusions in CZT and the role of photon recycling in TIBr and CZT. 
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I. INTRODUCTION 


A. NUCLEAR SECURITY THREATS 

The United States faces significant threats in the area of nuclear security. The 
threats come from both state and non-state actors, each of whom may employ different 
methods of attack and delivery systems resulting in differing levels of consequences. 

Potential state sources of nuclear weapons or fissile material include Iran, Russia, 
North Korea, and Pakistan [1], Iran is currently a nation of key concern. A November 
2011 report by the International Atomic Energy Agency (IAEA) found “credible” 
evidence that “Iran has carried out activities relevant to the development of a nuclear 
device” [2], Iran, however, has continued to defy UN Security Council Resolution 1969 
ordering it to suspend uranium enrichment operations [3], [4]. Terrorist groups have also 
actively sought nuclear weapons and fissile material for use against the United States. In 
1994, for example, al Qaeda failed in an attempt to purchase uranium, but a 2004 report 
by the National Commission on Terrorist Attacks on the United States asserts “al Qaeda 
continues to pursue its strategic objective of obtaining a nuclear weapon . . . [and] 
remains interested in using a radiological dispersal device” [5], 

Nuclear terrorism by a rogue state or group could take several forms. A group 
could employ a nuclear weapon or threaten its use in order to gain leverage. It may 
create an improvised nuclear device (IND) from illegally-obtained nuclear materials. It 
might release a radiological dispersal device (RDD) from radioisotopes designed for 
medical, research, or commercial use. Additionally, a terrorist cell could target 
installations that use radioactive materials, such as nuclear power plants [6]. 

Depending on the type of nuclear attack, the outcome could be catastrophic. One 
threat, with low probability but high consequence, is the detonation of a nuclear weapon 
in a major U.S. city. A ten-kiloton weapon detonated at Grand Central Station in 
Manhattan, New York, could potentially result in 500,000 instant fatalities, with 
hundreds of thousands of additional injuries and economic loss of over $1 trillion [1], 
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In order to launch an attack against the United States using a nuclear weapon, 
IND, or RDD, a terrorist group may need to smuggle the weapon into the U.S., or 
smuggle in the weapon components and then assemble it in the U.S. The weapon or its 
components could be hidden in a vehicle, train, or ship, concealed in personal luggage at 
an airport, or even walked across an area of unprotected border [7], [8], 

B. DETECTION CHALLENGES 

There are many complex and intertwined challenges to detecting illicit fissile or 
radioactive materials to protect against the nuclear threat. A few of these challenges are 
considered here. 

The vast borders of the United States and the sheer volume of cargo entering the 

U.S. each year confound radiation detection. In 2004, Admiral James Loy, former 

Deputy Secretary of the Department of Homeland Security, summarized the formidable 

scope of the problem facing the United States [9]: 

We must secure nearly 7,500 miles of land border with Canada and 
Mexico, across which more than 500 million people, 130 million motor 
vehicles, and 2.5 million rail cars pass every year. We also patrol almost 
95,000 miles of shoreline and navigable waters, and 361 ports that see 
8,000 foreign flag vessels, 9 million containers of cargo, and nearly 200 
million cruise and ferry passengers every year. We have some 422 
primary airports and another 124 commercial service airports that see 
30,000 flights and 1.8 million passengers every day. There are 
approximately 110,000 miles of highway and 220,000 miles of rail track 
that cut across our nation, and 590,000 bridges dotting America’s biggest 
cities and small towns. 

As of 2002, U.S. Customs was only physically inspecting 2%of the cargo entering the 
United States each year [10], Due to increased awareness of the nuclear threat, by 2006 
approximately 80% of the incoming cargo at port facilities was screened [11]. 

Detection and identification of radioactive materials depend on many factors, 
including the material quantity and decay scheme, the detection method (passive 
detection or active interrogation), the time allowed for detection, the distance between the 
source and the detector, the presence of shielding, and the knowledge of the operator. In 
passive detection, the detector simply detects radiation emitted from a source. In active 
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interrogation, however, a neutron or photon source bombards the source material, and the 
detector detects the secondary radiation that the source emits as a result of the 
bombardment [12]. Highly enriched uranium (HEU) is one of the materials of greatest 
concern. Unfortunately, it is also one of the most difficult to detect by passive means. Its 
low level of radioactivity and the nature of its radiation make it relatively easy to shield 
and therefore avoid detection [8]; depending on the purity of the HEU and the presence 
of shielding, hours of monitoring may be required to detect the HEU. In contrast, there 
are many legitimate sources of radiation that are easily detectable through passive means, 
including radiopharmaceuticals and many commercial goods such as ceramic tiles, 
fertilizers, and kitty litter. These legitimate sources can trigger radiation monitors, 
causing false alanns and costly delays [10], [13], [14]. 

In screening for radiation, the U.S. must balance two competing requirements - 
high sensitivity and high throughput. High sensitivity is required to distinguish between 
illicit radioactive materials such as HEU and legitimate sources such as 
radiophannaceuticals. High throughput is required so as not to impede the flow of 
commerce [13]. The Government Accounting Office (GAO) has documented cases 
where officials at U.S. ports of entry have sacrificed sensitivity for throughput. In one 
instance, the sensitivity of portal monitors was reduced in order to limit the frequency of 
nuisance alarms from naturally occurring radioactive materials. In another case, officials 
allowed trucks to pass through portal monitors at speeds greater than optimal for the 
detection of nuclear material [8], [13]. These cases highlight the magnitude of the 
detection problem. By 2006, the U.S. had screened 80 million cargo containers with 
radiation portal monitors, which triggered 318,000 false alarms. Most of the false alarms 
came from naturally occurring radioactive materials, and none of the alarms were 
attributed to illicit radioactive materials. Advanced detectors are therefore required that 
will allow for more accurate detection of illicit radioactive materials with fewer false 
alarms [15]. 

Current detection technology fails to meet the new demands of homeland defense. 
National security applications require ruggedized instruments that can rapidly measure 
the energy and identify the type of incident radiation from gamma and neutron sources. 
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For most applications, detection of gamma rays with energies between 10 keV and 
10 MeV is needed. In the lab and some limited field applications, cryogenically cooled 
high purity germanium (HPGe) semiconductor detectors can identify the isotopic 
signature of gamma rays with the required sensitivity and spectral resolution. The 
cooling requirement, however, excludes the use of FIPGe for many national security 
needs. Instead, many hand-held radionuclide identification devices use inorganic 
scintillators such as thallium-doped sodium iodide and many radiation portal monitors 
used at borders to scan vehicles use scintillating plastics such as polyvinyl toluene. 
These HPGe substitutes have poor energy resolution, as shown in Figure 1, which makes 
identification of detected radionuclides difficult. Neutron detection for neutrons arising 
from fissile material is also required. Currently, however, there is no commercially 
available portable instrument for rapidly identifying neutron spectra. Advanced radiation 
detector materials are therefore needed to meet the needs of national security [9], [16]. 



Figure 1. Comparison of energy resolution for natural background using high purity 

germanium, thallium-doped sodium iodide, and plastic scintillators. After [13]. 

C. MATERIALS PROBLEM 

Two primary classes of solid-state radiation detection materials are scintillators 
and semiconductors. With both inorganic scintillators and semiconductors, incident 
radiation produces electron-hole pairs (EHPs) in the materials. For scintillation detectors, 


4 












these EHPs are further converted to optical photons. These photons traverse the material 
to the electrodes, where they are collected and converted to an electrical signal. For 
semiconductors, the EHPs themselves are the final carriers of information; under the 
application of a strong electric field, these free carriers move to the electrodes, where 
they produce a current pulse that is proportional to the number of charge carriers 
collected. In general, semiconductor detectors have better resolution than scintillation 
detectors. This is due in part to the direct conversion of EHPs to electrical signal in 
semiconductors, without the intennediary step of converting EHPs to photons as required 
in scintillation detectors [16]. 

Decades of research and development are required to grow crystals with suitable 
properties for use as semiconductor radiation detectors. Researchers have investigated 
materials that can deliver the same energy resolution as Ge while operating at room 
temperature as far back as the introduction of lithium-drifted Ge detectors in the 
1960s [17]. Typically, however, new semiconductor materials have only been discovered 
at a rate of approximately once per decade [16]. This pace is inadequate for the needs of 
national security. 

One of the challenges in characterizing and developing new radiation detection 
materials is that key indicators of detector performance, such as charge carrier lifetime, 
are not intrinsic properties of the material. Instead, carrier lifetime depends on the nature 
of both the impurities and native defects in the material, which in turn depends on the 
crystal growth process. First principles theories therefore cannot be used to reliably 
predict suitable detector material candidates. Instead, each candidate material must 
undergo a painstaking experimental process to develop a viable crystal growth method, 
improve the material properties to allow measurement and characterization, and refine the 
material properties to ensure high perfonnance detector operation [17]. 

In the past, the pace and effectiveness of discovery of new candidate materials for 
radiation detection have been limited by a tendency to validate material performance only 
through fabrication and testing of actual radiation detectors [16]. New methods are 
therefore necessary that can probe the material properties of interest without requiring the 
cost and time burden of full device fabrication. 
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D. RESEARCH METHOD AND GOALS 


To accelerate the development of existing materials and the discovery of new 
materials for gamma ray detectors, there are several critical research needs. Two of these 
needs are improved characterization of detector material properties and improved models 
of radiation-detection physics in detector materials [16]. This research aims to contribute 
to both. 

The critical national security need for gamma ray detection is the ability to 
identify individual radionuclides. Each radionuclide emits a unique distribution of 
gamma radiation energies. To identify a radionuclide, high resolution energy 
spectroscopy is therefore essential [16]. For high energy resolution, a semiconductor 
detector needs uniform charge transport. A key measure of charge transport is the 
mobility-lifetime (/lit) product, which is discussed in greater detail in Chapter II. A 
uniform //r product across a detector results in high energy resolution, whereas variations 
in the //rproduct degrade energy resolution [18]. 

This research focuses on quantifying variations in the //r product in bulk 
semiconductors through materials characterization and modeling. A micro¬ 
characterization technique called transport imaging is being developed at the Naval 
Postgraduate School. Transport imaging is a contact-free all-optical technique that 
allows one to determine the //r product of a sample, without the need for contact 
fabrication or device processing. 

The broad objective of this research is to use a 2D carrier diffusion model in TIBr 
and CZT to demonstrate the application of transport imaging as a micro-analysis tool for 
characterizing the charge transport properties in bulk semiconductors for nuclear 
radiation detection. The specific goals of this research are to: 

• Advance the development of a rapid assessment micro-characterization 
technique to determine the magnitude and spatial variation in the /lit 
product at 2 pm resolution. 

• Investigate the dependence of the //r product on temperature and doping in 
TIBr and CZT. 


6 



• Contribute to the development of a defect level map of Se and Pb in TIBr. 

• Investigate the impact of Te inclusions on charge transport in CZT. 

• Investigate the effects of photon recycling on diffusion length in TIBr and 
CZT. 
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H. M. Smith, D. J. Phillips, I. D. Sharp, J. W. Beeman, D. C. Chrzan, N. M. 
Haegel, E. E. Haller, G. Ciampi, H. Kim, and K. S. Shah, “Electronic effects of Se 
and Pb dopants in TIBr,” Appl. Phys. Lett., vol. 100, no. 20, pp. 202102-1 - 
202102-4,2012. 


K. E. Blaine, D. J. Phillips, C. L. Frenzen, C. Scandrett, and N. M. Haegel, 
“Three-dimensional transport imaging for the spatially resolved detennination of 
carrier diffusion length in bulk materials,” Rev. Sci. Instrum., vol. 83, pp. 043702- 
1 -043702-7,2012. 
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“Spatial variation in mobility-lifetime product in bulk TIBr and CZT,” Presented 
at APS March Meeting, Boston, MA, 1 March 2012. 
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II. RADIATION SPECTROSCOPY AND THE jut PRODUCT 


A. GAMMA-RAY RADIATION SPECTROSCOPY 

When a gamma-ray enters a semiconductor radiation detector, an “energy 
cascade” is initiated, converting a single quantum of radiation through a series of 
interactions with electrons and atomic nuclei to a distribution of electrons and holes that 
produce a current signal in an instrument readout system. The cascade begins when an 
absorbed gamma-ray creates primary electrons, predominantly through photoelectric, 
Compton, and pair production interactions. These primary electrons then undergo a 
series of quantum energy transfer processes, producing electrons, photons, plasmons, and 
phonons [16]. The cascade ultimately results in thermalized EHPs with energies on the 
order of 1 eV, phonons, and atoms displaced by elastic energy transfers. Under the 
application of a strong electric field, these thennalized EHPs are swept out of the detector 
medium and recorded as a current pulse. The current pulse is proportional to the number 
of charge carriers collected, and therefore the energy of the initial radiation [16]. This 
process is the basis for gamma-ray spectroscopy. 

If the charge collection process is incomplete due to limited carrier transport 
and/or carrier trapping, a lower current pulse will be recorded, resulting in an inaccurate 
evaluation of the energy of the initial radiation. This translates into a broadening of the 
peak in the energy spectrum and lower detector resolution [19]. 

B. MOBILITY-LIFETIME PRODUCT 

The jut product is a key measure of charge transport in a semiconductor. It is 
comprised of two components—the drift mobility, //,/, (commonly referred to as simply 
the mobility, //), and the recombination lifetime, r, (commonly referred to as simply the 
lifetime, z). The mobility determines the drift velocity, vv, which is the speed at which 
charges move through a material under the application of an electric field, E, through the 
relationship 
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v d =fiE 


(1) 

The lifetime is the average time that carriers exist before recombining, becoming trapped, 
or scattering. Together, the /ur product determines the drift length, X, and diffusion 
length, L, for free carriers in a semiconductor. The drift length is a measure of the 
average distance that free carriers move though a material due to an applied electric field 
before recombining. The drift length is given by 

X- /lit E (2) 

The diffusion length is a measure of the average distance that charges move through a 
material due to diffusion before recombining. The diffusion length is given by 

where ks is the Boltzmann constant, T is the temperature in K, and e is the fundamental 
unit of charge. 

C. DETECTOR RESOLUTION AND CHARGE COLLECTION 


The resolution of a gamma-ray detector is commonly defined as the energy peak 
full-width at half maximum (FWHM, or AE) divided by the centroid of the energy peak 
E 0 , given by 


R = 


FWHM _ AE 
E, F 0 


( 4 ) 


By convention, resolution is quoted as a percent [20]. An idealized detector response and 
FWHM is shown graphically in Figure 2. 
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Figure 2. Simulated gamma-ray peak with FWHM and centroid energy E 0 . After [20]. 
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The resolution of a detector depends on the FWHM, and the FWHM depends 
primarily on three factors—carrier statistics, electronic noise, and charge carrier 
collection. Each factor contributes to peak broadening in the energy spectrum. The 
FWHM is given by 

A E +{t£ s ) ! +(6E c f . (5) 

The first term, AEp, represents Fano noise and accounts for the inherent statistical 
fluctuations in the number of charge carriers created. The second term, AE e , represents 
electronic noise and is dominated by leakage or dark current. The final term, AE C , arises 
from incomplete charge collection [20]. 

Figure 3 shows the relative contribution to peak broadening from each term for an 
HPGe detector with a volume of 86 cm . The greatest contribution comes from 
incomplete charge collection [20], Variations in the ju t product within the detector affect 
this term. 



Figure 3. Contribution to FWHM due to Fano noise (AE e ), 
electronic noise (AE e ), and charge collection ( AE C ). 
After [20], 
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The influence of variations in the /uz product on detector resolution can be seen 
indirectly by considering the charge collection efficiency (CCE). The CCE is a ratio of 
the initial amount of charge Qo created by EHPs due to the absorption of a gamma-ray to 
the total charge induced Q on the electrodes of the detector, or 

CCE = Q~. (6) 

Qo 

Ideally, a detector would have a CCE of unity, meaning all created charge is 
collected. In real semiconductors, however, carrier diffusion and carrier trapping result 
in the loss of carriers and reduced charge at the electrodes [21]. Statistical models and 
curve-fitting methods have been extensively developed to detennine the CCE [22], [23]. 

One common curve-fitting technique for determining the actual CCE is using the 
Hecht method. A radiation source such as a laser or an alpha emitter is used to generate 
charge near the detector cathode or anode. The induced charge signal is then measured as 
a function of the bias voltage across the detector [24], The /uz product is obtained by 
curve fitting to the Hecht equation, given by 


CCE = — J1 - exp ^ 

L , 2 . 


1 , ( x n 

+—< 1 - exp —- 

L . 2./, 


where L is the detector thickness, xq is the distance from the cathode to the point of 
charge creation, and A e and A/, are the drift length of electrons and holes, respectively. 
The drift lengths are 

^,=K £ (8) 


\ ={f* T ) h E ( 9 ) 

where (juz) e and (/uz) h are the /lit products of electrons and holes, respectively. Small 
ratios of AIL for either electrons or holes reduce the CCE, resulting in peak broadening in 
the energy spectrum [21]. Spatial variations in the /uz products also cause variations in 
these ratios, also contributing to peak broadening. 
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The influence of variations in the jut product on detector resolution can be seen 
directly in certain cases. One case is considered here. Typically, the contribution of 
charge collection to peak broadening has been determined by the semi-empirical formula 

A E = sj(AE F ) 2 +(AE E f +(2.355) 2 a x E ai (10) 

where ai and a? are semi-empirical constants determined by best-fitting [21]. Kozoresov 
(2005), however, has shown that for pixelated detectors, 



where the trap factor G(E) can be given in closed form as 

f l 0 (e)) 2 

G(E) = (12) 

\ ) 

where L 0 (E) is the photon absorption length [25]. Here, it is possible to see directly that 
variations in the jut product would result in variations in the trap factor, which contributes 
to peak broadening. 

D. MOBILITY-LIFETIME PRODUCT MEASUREMENT TECHNIQUES 

Several techniques are available for measuring the material unifonnity and 
variations in the jut product of a semiconductor. These methods include alpha particle 
response characterization [26], high resolution mapping techniques using synchrotron X- 
rays, nuclear microprobes and ion microbeams, [27] - [29] and transport imaging [30]. 
Except for transport imaging, all of the mapping techniques require device processing 
and some form of the Hecht equation to extract the jut product. 

One limiting assumption of the Hecht relation is that charge trapping is uniform 
within the bulk of the material and detrapping is negligible [19], [31], Real crystals, 
however, have inhomogeneities in charge trapping. Additionally, the jut products 
obtained using the Hecht relationship are average values over the intercontact distance. 
Therefore, although a map of the jut product is generated, the /lit products represent 
averages over a large distance, with distances ranging from millimeters to centimeters or 
more, depending on the type of detector. While this is relevant for detector performance, 
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it is limited in advancing the type of micro-scale characterization that will be required for 
more fundamental understanding of material properties and the actual role of various 
defects and dopants that limit performance. 

Transport imaging is a unique imaging method based on the use of a scanning 
electron microscope (SEM) that allows measurements of the diffusion length and 
extraction of the jlit product in luminescent materials without the need for device 
processing, applied electric fields, or averaging over large intercontact distances. In 
transport imaging, the SEM generates free charge carriers in a semiconductor, and a 
charged coupled device (CCD) array camera records the motion of these charge carriers 
as they diffuse by imaging the distribution of their recombination luminescence. The 
technique has been applied previously to obtain /jt products in bulk GaAs, minority 
carrier lifetimes in heteroepitaxial GalnP, and minority-carrier diffusion lengths in 
ordered GalnP [32] - [35]. This work will advance the application of transport imaging 
as a micro-analysis tool for the characterization of bulk materials for nuclear radiation 
detectors. 
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III. TRANSPORT IMAGING 


A. EXPERIMENTAL SETUP 

Cathodoluminescence (CL) is the emission of light from a material under 
excitation by an electron beam. In order to apply transport imaging to materials, the 
samples must have a CL signature. For CL, samples are placed on the continuous-flow 
cold stage of a JEOL 840A scanning electron microscope (SEM) and may be cooled with 
liquid helium. Temperatures from 300 K to 5 K can be obtained. The electron beam is 
incident normal to the sample. In this work, the SEM beam excitation energies varied 
between 10.0 keV and 20.0 keV, the probe currents ranged from 3><10' 10 A to 6 X 10' 10 A, 
and the magnifications ranged from 1000x to 2000*. An Oxford CL system with a 
0.25 m path length monochromator was used to perform spectroscopy. The detector is a 
thermoelectrically cooled photomultiplier with response from 300-900 nm. 
Spectroscopy can be performed for emission obtained with the SEM operating in either 
spot or picture (area) mode. It is also possible to acquire CL area maps in either 
panchromatic or monochromatic mode. The SEM and CL system are shown in Figure 4. 



Oxford CL System 


SEM Chamber 


CCD Camera 


SEM Electronics 


Figure 4. Experimental set up for CL and transport imaging using the electron 
beam from the SEM and an external CCD array camera. 


For transport imaging, an optical microscope is inserted into the vacuum chamber 


of the SEM using a retractable arm. A small hole, or objective pipe, in the microscope 
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allows passage of the electron beam. A 45 degree mirrored surface collects the 
luminescence and directs it through a light guide to a 1284 pixel x 1472 pixel 
thermoelectrically-cooled Si CCD array camera external to the SEM chamber. The pixel 
size on the CCD camera is 6.8 pm x 6.8 pm, and the magnification of the optical 
microscope is ~20 x , resulting in optical images with spatial resolution of -400 nm/pixel. 
The setup for transport imaging is shown in Figures 4 and 5. 



Figure 5. System schematic for transport imaging using the electron beam 
from the SEM and an external CCD array camera. 

B. MATHEMATICAL MODEL OVERVIEW 

In transport imaging, the electron beam of the SEM generates the EHP charge 
carriers in the semiconductor material. The charge carriers can diffuse through the 
material and recombine at different locations, emitting light. The CCD array camera 
records the location and intensity of the recombination luminescence. By fitting a model 
to the recorded intensity distribution, it is possible to extract the diffusion length and 
calculate the /lit product. 

The model (ID, 2D, or 3D) employed in transport imaging to extract the diffusion 
length depends on the semiconductor material geometry (nanowire, thin film, or bulk 
sample) and the SEM excitation mode (spot or line mode). In the following section, we 
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review the models for the study of nanowires and thin films, and then describe the 
extension to bulk samples that is applied in this work. The semiconductor material 
geometries, SEM excitation modes, and models are summarized in Figure 6. 



Figure 6. Summary of the semiconductor material geometries (nanowire, thin 

film, or bulk sample), SEM excitation modes (red dot = spot mode, red line = line 
mode), and models (ID, 2D, or 3D). Arrows indicate diffusion directions. 


For nanowires, the SEM beam is operated in spot mode, and a ID model is applied. 
The intensity of recorded light / is given by 

-f-1 

/~e Uj (13) 

where x is the horizontal distance from the SEM beam. The ID model for transport 
imaging has been applied previously to GaN nanowires using near field optical 
microscopy [36]. 
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For thin films, the SEM beam may be operated in either spot mode or line mode, 
providing either a quasi point source or ID line source for carrier excitation. For spot 
mode, a 2D diffusion model is applied, and the intensity of recorded light is given by 

f r ^ 

I~K 0 - (14) 

\ L J 

where Kn is a modified zeroth order Bessel function of the second kind, and r is the radial 

distance [r = sjx 2 + y 2 j from the point source in the plane of the thin film to an arbitrary 

location (x,y) where recombination occurs and light is generated. The 2D model for 
transport imaging in thin films has been applied previously to the study of ordered GalnP 
for solar cell applications [35]. To take advantage of symmetry, the SEM may also be 
operated in line mode for thin films. This reduces a 2D problem to a ID problem. For 
line excitation in thin films, a ID model is applied, and the intensity of the recorded light 
is given by Equation (13). The ID excitation model for transport imaging in thin films 
has been applied previously to the study of dislocation networks in heteroepitaxial GalnP 
[34], [35], 

For bulk samples, the SEM beam may be operated in either spot mode or line 
mode. For spot mode excitation, a 3D diffusion model is applied, and the intensity of 
recorded light is given by 



where r is the radial distance = ^x 2 + y 2 + z 2 j from the point source on the surface to 

an arbitrary location (x,y,z) in the bulk sample where light is generated [37]. By again 
taking advantage of symmetry, the problem may be reduced from 3D to 2D by operating 
the SEM in line mode. For line mode excitation, a 2D diffusion model is applied, and the 
intensity of the recorded light is obtained by integration through the bulk (z-direction) in 

Equation (14), with the modification that r = ^x 2 + y 2 + z 2 . The 2D model for transport 
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imaging in bulk samples has been applied for the first time to GaAs by Blaine [32], [33]. 
For this present work in TIBr and CZT, the 2D model is also applied, as described in 
greater detail below. 

C. BULK SAMPLE 2D MATHEMATICAL MODEL AND ASSUMPTIONS 

1. Mathematical Model 

A mathematical model has been developed by Blaine [32], [33] to extract the 
diffusion length L and the mobility-lifetime product //r from the luminescence 
distribution using the 2D model for transport imaging in bulk material. The model is 
briefly described here. 

In the 2D model for bulk samples, the SEM is used to generate a line of charge 
along the y-direction of the sample, as shown in Figure 7(a). The rate of charge 
generation per unit length is given the symbol g [charge/cm-s]. When electrons impact 
the bulk sample, they create a generation volume of EHP charge carriers [EHPs/cm ] 
below the surface of the sample, primarily through photoelectric, Compton, and pair 
production interactions. The centroid of the generation volume is modeled as a delta 
function at the position x = 0 and depth z = z 0 below the surface, as shown in Figure 7(b). 



Figure 7. (a) The SEM generates a line of excess charge carriers across the 

bulk sample. The bulk sample is modeled as a series of semi-infinite planes, 
labeled A through N. (b) For each plane, the line source is modeled as a delta 
function at depth zo below the surface. The carrier concentration u(x) is 
determined for each point on the x-axis of the plane, labeled x = ai ... on- 
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The charge carriers created in the generation volume diffuse and recombine. The 
concentration of charge carriers at any location r in the material is given by the symbol 
u ( f). In the 2D model, both the line source of excitation and the bulk sample are 

modeled as infinite in the y-direction. This is a valid assumption as long as the length of 
line source is much greater than the carrier diffusion length in the material. Under this 
condition, there is no change in the carrier concentration along the y-direction. As a 
result, 

w(r) = u(x,z). (16) 

With this simplification, the bulk sample can be treated as set of semi-infinite planes, 
labeled A through N in Figure 7(a), with diffusion in the x and z directions only. 
Analysis can therefore be restricted to a single plane. Analysis of the SEM line source on 
the bulk sample has therefore been reduced to analysis for a point source in a plane. 

The diffusion of charge carriers is governed by the continuity equation. The time 
rate of change of the charge carrier concentration depends on the rate of charge carrier 
generation by the SEM beam, the rate of recombination, and carrier diffusion, given by 

| ] -u(r,t) = G-R + -(W-j) (17) 

where u[r,t ) is the charge carrier concentration at position rand time t, G is the charge 
carrier concentration generation rate, R is the charge carrier concentration recombination 
rate, e is the fundamental unit of charge, and J is the current density. 

For transport imaging with the beam operating in line mode, the time required for 
the beam to raster back and forth is short compared to the selected exposure time, and the 
integrated dwell time for the beam at any point is very long compared to the carrier 
lifetime. This creates steady-state conditions for the charge carrier distribution. Under 
steady-state conditions, 

—u(r,t) = 0. (18) 

dt v ’ 

The charge carrier generation rate is given by 

G = g-S(x)S(z-z 0 ). (19) 
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(20) 


The charge carrier recombination rate under steady-state conditions is 



where z is the recombination lifetime. The current density in the absence of an electric 
field is 

J = eDVu(r ) (21) 

where D is the ambipolar diffusion constant. The ambipolar diffusion constant is related 
to the relevant diffusion length through the relationship 

L = 4Dt . (22) 

Combining Equations (16) through (22) results in 

V 2 w(x,z)--^n(x,z) + — S(x)S(z-z 0 ) = 0. (23) 

L D 

Two boundary conditions are specified. The first boundary condition forces the 
excess carrier concentration to zero at large distances from the source. This condition is 
given by 

u{x,z) —>0 as yjx 2 +z 2 —»oo. (24) 

The second boundary condition introduces the surface recombination velocity S and 
specifies that the diffusion of charge carriers normal to the surface is proportional to the 
charge carrier concentration near the surface. This condition is given by 


D— m(x,z) = Su[x,z^ atz = 0. 


In order to fit the model to the luminescence data collected through transport 
imaging, it is necessary to solve Equation (23) for u(x,z) subject to the boundary 
conditions given by Equations (24) and (25). After solving for u(x,z), it is necessary to 
integrate w(x,z) through the depth of the material, from z = 0 to infinity, since below band 
gap photons can escape from the bulk of the material. This results in a model solution for 
the excess carrier concentration u(x) for each point on the x-axis of the plane, labeled x = 
ai ... a n for plane A in Figure 7(b). The analytic model solution for the integrated excess 
carrier concentration, as derived by Blaine in [32], is 
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(26) 


f(x) = A 




JX (A t 2 - - dt+ Jk, 


yLj 


yjt 2 +x 2 


z o V 


L 




L 


dt 


where K 0 is the zeroth-order modified Bessel function of the second kind, Kj is the first 
order modified Bessel function of the second kind, and 

A = ^A. (27) 

n D 

The model solution in Equation (26) is based on the assumption that the recorded 
luminescence intensity is proportional to the integrated carrier concentration u(x). The 
model is fit to the normalized luminescence intensity data to obtain the diffusion length. 
The fitting is done using the Levenberg-Marquardt algorithm to minimize the square 
residual [32], [33], 

The diffusion length is related to the /ur product by combining Equation (22) with 
the Einstein relationship, given by 


D = ^ (28) 

e 

where ks is the Boltzmann constant and T is the temperature in K. Combining Equations 
(22) and (28) results in 


L = ■ ( 29 ) 

Equation (29) allows for the extraction of the //r product from the diffusion length 
obtained from the fit of Equation (26). 


2, Model Assumptions 


Two key assumptions of the model are low injection conditions and constant 
lifetime. When excess EHPs are created in a semiconductor, the spontaneous radiative 
lifetime for band to band recombination is given by 


t r- - 

R R 


NN 


Ak B T(n 0 +Po) 


-i-i 


1 + 


\% +Pi 


0 J 


(30) 
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where n e is the excess carrier concentration, R is the recombination rate, N c is the density 
of states at the conduction band edge, N v is the density of states at the valence band edge, 
A is the Einstein A coefficient, ks is the Boltzmann constant, T is the temperature, no is 
the equilibrium electron concentration, and po is the equilibrium hole concentration. 
Under low injection conditions in /7-type material, n e « n 0 , resulting in 

(31) 


( N N 

r i ] 

Uw 

Uo +P 0 ) 


Typically, the coefficients of the carrier concentrations are combined into one coefficient, 
where 

Ak a T 


B = 


NN 


(32) 


Combining Equations (31) and (32) yields 


1 


U? =- 


B (%+Po) ' 

For n-type doping, where no »po and no = N D , Equation (33) reduces to 

1 


(33) 




BN, 


(34) 


Thus, the model for detecting, for example, the minority carrier (hole) diffusion length in 
//-type material, is based on a single constant lifetime under low-injection conditions [38], 


D. 


TRANSPORT IMAGING EXAMPLE 


This section presents a detailed example of the transport imaging technique. A 
sample of TIBr doped with copper, iron, and zinc is used for demonstration. 

For transport imaging in a bulk sample, the SEM beam is operated in line mode to 
generate a line of charge carriers along the v-direction of the sample, as shown in the 
image of the recombination luminescence in Figure 8(a). To avoid edge effects, in this 
example a 10 pm segment is extracted from the center of the ~ 60 pm line scan image; 
segments of other lengths may be chosen, with 40 pm segments most commonly selected. 
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Planar and oblique views of the recombination luminescence intensity distributions from 
the extracted 10 jam segment are shown in Figures 8(b) and 8(c). Note that since the 
resolution of the imaging system is 0.4 pm/pixel, the 10 pm segment is comprised of 
25 individual pixels along the y-di recti on. This results in 25 separate intensity 
distributions in the 10 pm segment. 

The diffusion length information is contained in the intensity distribution 
perpendicular to the line scan, along the x-axis. Each of the 25 distributions shown in 
Figure 8(c) or 8(d) is normalized, and then averaged to create Figure 8(d). Often, there is 
a small bump to the left of the peak due to a reflection in the optical system, although this 
bump is not visible in visible in Figure 8(d). Therefore, by convention, data from only 
one side of the peak (x > 0) in Figure 8(d) are used for least squares fitting to the model. 

Figure 8(e) shows the normalized mean intensity distribution data and the best-fit 
line from the model given by Equation (26). For this example fit, data were taken 
beginning at x = 4 pm to avoid effects of the generation volume. The model provides the 
diffusion length as one of the fitting parameters. This diffusion length is then used in 
Equation (29) to calculate the //rproduct. 
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Figure 8. (a) SEM line scan luminescence image of TIBr crystal. 

(b) Planar and (c) Oblique views of a 10 jam segment (25 pixels) of the 
recombination luminescence intensity distribution. The diffusion length 
information is contained in the intensity distribution perpendicular to the line 
scan, along the x-axis. (d) Normalized mean recombination luminescence 
intensity distribution for the TIBr crystal, (e) Data from right side of distribution 
from Figure 4(d) and best-fit line from least squares fitting to the model. 
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IV. INITIAL INVESTIGATION OF BULK THALLIUM BROMIDE 


A. OVERVIEW OF THALLIUM BROMIDE 
1. Properties of TIBr 

TIBr is an attractive compound semiconductor alternative to HPGe for room 
temperature radiation detection and gamma spectroscopy [39]. Its high density 
(7.56 g/cm ) and atomic number (Zti = 81, Zb,- = 35) result in high photon attenuation and 
detection efficiency [16], [40]. The wide indirect bandgap (2.68 eV) and high room- 
temperature resistivity (~1() 10 Q cm) of TIBr reduce leakage current and shot noise [41], 
[21]. TIBr has a simple cubic crystal structure and a low melting point (460 °C), which 
makes it easy to grow single crystals [42], Advances in TIBr purification techniques and 
growth processes have led to steady improvements in the mobility-lifetime products of 

3 2 

electrons and holes, ( jur) e and (jUT) h , respectively. Values of ( jur) e of 6.5 X 10' cm /V and 

{jur)h of ~4x 10" cm /V at 300 K have recently been reported that are comparable to those 
of cadmium telluride (CdTe) and Cadmium Zinc Telluride (CZT). With cooling to 18° C 
and depth correction techniques, energy resolutions of < 1% (FWHM) for the 662 keV 
peak of 137 Cs have been achieved for TIBr array detectors [40]. Table 1 summarizes 
select properties of TIBr. 
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Table 1. Properties of TIBr. 


Property 

Value 

Reference 

Atomic Number 

Z T i: 81, Z Br : 35 

[43] 

Density 

7.56 g/cm 3 

[43] 

Bandgap 

2.68 eV 

[43] 

Resistivity 

~10 10 tlcm 

[43] 

Electron Mobility 

40 (cm 2 /V-s) 

[44] 

Hole Mobility 

12 (cm 2 /V-s) 

[44] 

Electron /ur Product 

6.5xlO J cm 2 /V 

[40] 

Hole /ut Product 

~4xl0" 4 cnUV 

[40] 

Melting Point 

460 °C 

[21] 

Crystal Structure 

Cubic (CsCl) 

[21] 


Despite recent advances in crystal growth and purification, challenges remain 
before TIBr can be widely deployed as a room temperature nuclear radiation detector. It 
displays significant ionic conductivity under bias, which degrades detector performance 
over time. Additionally, the origin of its good transport properties ( jut product) and high 
resistivity are not well understood [45], [46]. 

Semiconductors such as TIBr with high resistivity (> ~10 7 Q-cm) are classified as 
semi-insulating [47]. High resistivity is generally achieved through the intrinsic nature of 
the material, purification, or compensation doping. Most materials are far from intrinsic, 
and purification is required. Ge has been made semi-insulating through purification of 
electrically active dopants to the level of ~ 10 cm' and by cooling. For most materials, 
however, such pure purification levels are not achievable, and compensation of existing 
impurities is necessary. Dopants may be introduced or defects may be controlled in order 
to create deep levels that compensate shallow donors or acceptors and “pin” the Fermi 
level near the center of the bandgap [48], [49]. Unfortunately, the same deep levels that 
ensure high resistivity are also often highly efficient recombination centers and carrier 
traps, which can reduce the //r product of the material [44], [50], 

In TIBr, the origin of the high resistivity is not well understood, and competing 
models have been proposed. Du (2010) has shown theoretically that low energy native 
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defects do not induce deep electron traps in TIBr, and that the dominant native defects in 
TIBr are Schottky defects (pairs of T1 and Br vacancies) which pin the Fermi level near 
the middle of the bandgap [45]. Elshazly (2010), however, argues that the pinning of the 
Fenni level in TIBr is due to the presence of deep level traps that compensate shallow 
level traps [44]. The contributions of ionic conductivity must also be fully investigated in 
order to understand and isolate the electronic contributions. 

TIBr is also not a well-characterized semiconductor material, and much more 
work remains to be done to develop high volume, well-controlled growth and 
performance. Because transport imaging isolates only electronic transport, it can play a 
role in this work. 

2. Transport Imaging in High Resistivity TIBr 

As a semi-insulating material with the Fermi level pinned near the middle of the 
gap, the concentration of electrons and holes in TIBr is very low and approximately equal 
[32], [43], The luminescence that we measure in transport imaging will therefore be due 
to recombination of the EHPs that we created with the SEM beam. In transport imaging, 
no electric field is applied to the sample. The EHPs therefore diffuse together as 
excitons. The diffusion of charge carriers is limited by the mobility of the slower carrier, 
which is the hole in TIBr. The diffusion length measured is therefore, to first order, an 
ambipolar diffusion length. The mathematical fit to the model subsequently generates an 
ambipolar //r product. 

B. CRYSTAL GROWTH 

Analysis in TIBr began with sample 45-S10 from crystal ZR26-45YD. The 
detector-grade crystal was grown by Radiation Monitoring Devices, Inc. (RMD) for 
making planar and array detectors. For crystal growth, a synthesized 5N purity TIBr 
compound was purchased from a commercial vendor and used as the starting material. 
The crystal went through 100 passes of zone refining and purification before being grown 
by the travelling molted zone method. It was grown in an open boat under a HBr and Ar 
mix gas flow. 
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c. 


GLOW DISCHARGE MASS SPECTROMETRY 


Glow discharge mass spectrometry (GDMS) was performed by Evans Analytical 
Group (EAG) to determine impurity concentrations. GDMS is an important 
characterization technique for TIBr because it can detect a wide range of impurities. The 
current state of TIBr is such that it can contain many different types of impurities in the 
concentration range of 10 15 — 10 19 cm' 3 . GDMS results for sample 45S-10 are shown in 
Figure 9. The impurities that remain in the crystal in the highest concentrations after 
zone refining and growth are aluminum (65 ppm wt, 1.1x10 cm' ), silicon (2.3 ppm wt, 
3.7x 10 17 cm' 3 ), and barium (3.3 ppm wt, 1.1x10 17 cm' 3 ). These residual levels indicate 
that sample 45-S10 is unintentionally highly doped. 
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Figure 9. GDMS results for TIBr sample 45-S10. 
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D. 


CATHODOLUMINESCENCE 


In order to perform transport imaging on a sample, the sample must have a CL 
signature. CL mapping and spectroscopy were performed on TIBr sample 45-S10 to 
determine if the samples would luminesce. Since TIBr is an indirect band gap material, it 
was not kn own if TIBr would emit sufficient luminescence. To our knowledge, only one 
previous report exists for CL in TIBr, and the spectra was obtained at 95 K and showed 
limited signal to noise ratio [51]. 

Room temperature panchromatic CL mapping was performed on sample 45-S10, 
as shown in Figure 10. The SEM operating parameters were 20.0 kV, 6 X 10" 10 A, and 
2000 x magnification. Due to the depth of penetration of the SEM beam, effects 
of surface contamination are negligible. Intensity variations are observed on a scale of 
~10 pm. It is this type of spatial variation that can be quantitatively studied with 
transport imaging. 



Intensity (Arb. units) 


Figure 10. Cathodoluminescence map (false color) of TIBr sample 45-S10. Intensity 

variations are observed on a scale of -TO pm. 
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CL spectra were measured for TIBr sample 45-S10 at 11 K and 295 K. For the 
CL system, the slits were set to 1 mm, and the step size was 1 nm. 

The CL spectrum at 11 K is shown in Figure 11. Four peaks, labeled A through 
D, are observed: A at 411 nm (3.02 eV), B at 461 nm (2.69 eV), C at 477 nm (2.60 eV), 
and D at 691 nm (1.80 eV). Peaks A and B are consistent with reported values for the 
direct and indirect gap of TIBr at 3.02 eV and 2.68 eV, respectively [52]. The origin of 
peak C at 2.60 eV is unknown. Peak C may represent a phonon replica of the LO phonon 
associated with peak B. The energy separation between peaks A and B is 90 meV. 
Inelastic neutron scattering experiments in TIBr have shown energy separations between 
phonon replicas of 13.9 meV. Shimizu (1974), however, has observed 26 meV 
separations between successive photoluminescence peaks in TIBr [53]. Peak C may 
therefore represent the third phonon replica of peak A. The spectrum also shows a broad 
peak at 691 nm (1.80 eV), which can be attributed to defect-related luminescence. 



Figure 11. CL spectrum for TIBr sample 45-S10 at 11 K. 


The CL spectrum at 295 K is shown in Figure 12. Gaussian fits to two peaks are 
included in the figure. Two features are noticeable in the room temperature spectrum. 
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First, a broad peak at 416 nm (2.98 eV) is now present in place of the sharp peak at 
461 nm (2.69 eV) seen in the 11 K spectrum. Second, the defect luminescence peak is 
now located at ~ 575 nm (2.16 eV). 

CL has been demonstrated in TIBr at both warm and cold temperatures. The CL 
peaks in TIBr are an unexplored area of research. To accurately assign all peaks, more 
detailed spectroscopy studies are required. For the purposes of transport imaging, 
however, there is a distinct luminescence signature, clearly attributable to bulk TIBr, 
which can be used to image carrier diffusion and recombination. 



Figure 12. CL spectrum for TIBr sample 45-S10 at 295 K. 


E. PHOTOLUMINESCENCE 

Photoluminescence (PL) was also performed on sample 45-S10. The PL system, 
located at the University of California, Berkeley, is comprised of a HeCd laser with an 
excitation wavelength of 325 nm, a UV/VIS spectrometer with a Si PMT detector, a 
385 nm cutoff filter, and a thennoelectric cooler. The laser was operated at maximum 
power of 50 mW, and measurements were recorded at a sample temperature of 10 K. 
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The PL spectrum and Gaussian fits to the peaks are shown in Figure 13. Six peaks, 
labeled A through F, were observed: A at 414 nm (3.00 eV), B at 428 nm (2.90 eV), C at 
475 nm (2.61 eV), D at 487 nm (2.55 eV), E at 503 nm (2.47 eV), and F at 644 nm 
(1.93 eV). 
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Figure 13. PL spectrum at 10 K for TIBr sample 45-S10. The dashed lines 

indicated Gaussian fits to the peaks. 


Peaks A and C are consistent with the direct and indirect bandgaps for TIBr. The 
origins of peaks B, D, and E are un kn own. The energy separation between peaks A and 
B is 98.0 meV. As suggested for the CL spectra, Peak B may represent a phonon replica 
of the LO phonon associated with peak A. Similary, peaks D and E are separated from 
peak C by energies of 64.3 meV and 145 meV, respectively. These peaks could therefore 
represent phonon replicas of peak C. It is also possible, however, that peaks B, D, and E, 
are not phonon replicas. They could instead be associated with defect-mediated 
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recombination. Peak F is associated with broad defect luminescence. Further doping 
experimentation and detailed spectroscopy are required to determine what levels the 
impurities introduce into TIBr. 

Figure 14 compares the CL and PL data for TIBr sample 45-S10. The CL and PL 
peaks at the direct gap (~ 3.0 eV) are closely matched. At the indirect gap, CL shows a 
peak at 2.69 eV (CL peak B), whereas PL shows a peak at 2.61 eV (PL peak C). The 
other peaks are more difficult to match. High resolution variable temperature CL and PL 
comparison studies are a prime area for future research. 
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Figure 14. CL and PL spectra of TIBr sample 45-S10 
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F. 


ENERGY DISPERSIVE X-RAY SPECTROSCOPY 


CL mapping shows spatial variations in the luminescence of the TIBr samples, but 
the technique does not give information as to any chemical origin of the variations. 
Energy dispersive X-ray spectroscopy (EDX) was performed to determine if the 
variations reflect non-unifonn distribution of the impurities. In EDX analysis, an 
electron beam bombards a sample, and a detector records the characteristic X-rays 
emitted in order to determine the elemental composition of the sample. EDX was 
performed at 15.0 kV, 2000 x magnification, and 295 K. 

EDX mapping of TIBr sample 45-S10 was performed. Figure 15(a) shows the 
distribution ofMai X-rays from thallium. Figure 15(b) shows the distribution of L a 1,2 
X-rays from bromine. Both distributions are uniform. This indicates that the thallium 
and bromine are unifonnly mixed, as expected, or that the non-uniformity is below the 
detection limit for EDX (0.1-1 at. %) [54]. 




1 Thallium: I QHi im ' Bromine: 

M (2.2683 keV) oUjJm L ttl 2 (1.4805 keV) 

Figure 15. EDX map for (a) thallium and (b) bromine. 

The characteristic X-ray spectrum for the TIBr sample is shown in Figure 16. 
EDX spectroscopy primarily detected characteristic X-rays from thallium and bromine. 
No other elements were detected above the EDX detection limits. This is consistent with 

19 3 

the GDMS results, which shows dopant/impurity levels below 10 cm". 
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Figure 16. EDX spectrum for TIBr sample 45-S10. 


G. TIME-OF FLIGHT SECONDARY ION MASS SPECTROMETRY 

In order to better understand the origin of the spatial variation observed in the CL 
maps, time-of-flight secondary ion mass spectrometry (TOF-SIMS) mapping was 
performed by EAG. The author traveled to the laboratories of EAG in Sunnyvale, CA, to 
observe the analysis. 

TOF-SIMS focuses a pulsed beam of Ga ions onto the surface of the sample. 
Secondary ions are generated by a sputtering process. Mass spectrometry then provides 
information about the molecular and elemental species present on the surface. The 
detection limit is 10 7 —10 10 at/cm 2 sub-monolayer. For mapping, the ion beam is rastered 
across the surface of the sample. The lateral resolution is -0.20 pm [55]. The samples 
were sputtered for 70 ns, with surface removal rate of-1 angstrom/second, resulting in a 
total material removal to a depth of -7 nm. Typically, the sputtering process removes 
contaminants on the surface. 

Figure 17 shows TOF-SIMS maps for thallium and bromine in TIBr sample 
45-S10. Figures 17(a) and 17(c) show maps for thallium and bromine prior to sputtering, 
and Figures 17(b) and (d) shown maps for thallium and bromine after sputtering. In each 
case, the maps of the unsputtered surfaces depict large intensity variations. This indicates 
that either the surface has contaminants that prevent detection of an underlying uniform 
thallium bromide matrix, or that the thallium bromide matrix itself is not homogenous, 
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with regions of high thallium or bromine concentration in one location, and regions of 
low concentration in another location. The maps of the sputtered surfaces are nearly 
uniform in their intensities, with some darker lines that may reflect scratches on the 
surface. The uniformity of the intensity indicates that the thallium bromide matrix is 
largely homogenous, and that the sputtering process has removed surface contamination. 




(d) 


% ):'■ • ■ '% 



%■ 


Br (sputtered) 10 \xm 


Figure 17. TOF-SIMS mapping showing elemental distribution of (a) thallium prior to 
sputtering, (b) thallium after sputtering, (c) bromine prior to sputtering, and 

(d) bromine after sputtering. 


TOF-SIMS maps for Na, Ca, K, Cu, F, and Cl are shown in Figure 18. For each 
element, an unsputtered and sputtered image is shown. Sputtered images for Na, Ca, and 
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K shown in Figure 18(b), (d), and (f) all depict intensity variations on a scale of 
~ 10-20 pm. This indicates nonuniformity in the distribution of these impurities. The 
intensity for Cu in Figure 18(g) and (h) is low, but the sputtered image appears to depict 
some nonunifonnity. For F and Cl, and sputtered images in Figure 18(j) and (i) also 
show low intensity variation, again indicating nonuniformity in impurity distribution. 


(a) 
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Na 
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Ca 
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Figure 18. TOF-SIMS mapping showing elemental distributions Na, Ca, K, Cu, F, and Cl 

before and after sputtering. 


H. TRANSPORT IMAGING—TEMPERATURE DEPENDENCE 

For transport imaging on TIBr sample 45-S10, the effect of SEM operating 
current first was investigated in order to determine the optimum operating parameters. 
Figure 19 shows a comparison of the intensity distributions at 20.0 keV for currents of 
lx 10' 10 A, 3x 10‘ 10 A, and 6 X 10“ 10 A for 100 lines averaged from a line scan. The 
distributions show good reproducibility for currents of 3 X 10' 10 A, and 6 X 10‘ 10 A. An 
operating current of lx 10' 10 A shows higher variations in the distribution beginning at 
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about 30 pm from the peak, possibly attributable to increased noise and background 
effects. Higher beam currents introduce undesirable surface effects due to local heating, 
whereas lower beam currents reduce the intensity to such a level that focusing the 
instrument is difficult. As a result, an operating current of 3><10' 10 A was found to be 
optimal, and it was used whenever possible. 



Figure 19. Evaluation of effects of beam current on intensity distribution 

at 20.0 keV in TIBr sample 45-S10. 


The effects of temperature on transport properties in TIBr sample 45-S10 were 
also studied. Figure 20 shows the normalized mean intensity distributions for 
10 different temperatures. The SEM operating parameters were 20.0 keV and 3 x 10‘ 10 A. 
The variation with temperature indicates that bulk material properties are playing a 
measurable role in the luminescence distribution. 
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Figure 20. Normalized mean intensity distributions for TIBr sample 45-S10 

for temperatures between 8 K and 102 K. 


Figure 21 shows a representative fit of the data at a temperature of 8 K. As 
highlighted in the semilog plot in Figure 21(b), the fit deviates from the data at ~ 15 pm 
from the source. The cause of this deviation has not yet been resolved, but it may be due 
to background subtraction limitations. 
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Figure 21. Fit to TIBr sample 45-S10 at 8 K shown on (a) linear and 

(b) logarithmic scale. 

Figure 22 shows plots of temperature vs. diffusion length and temperature vs. juz 
product obtained by fitting the data in Figure 20 to the model. The juz product decreases 
with increasing temperature over the range of 8 K-102 K. Owens (2003) and Onodera 
(2007) have reported on the temperature dependence of the /jz products for electrons and 
holes over temperature ranges of 265 K-278 K and 233 K-313 K, respectively [56], [57]. 
To the author’s knowledge, this work presents the first report of the ambipolar /vr product 
over the range of 8 K-102 K. 
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Figure 22. Temperature dependence of (a) diffusion length and (b) //r product for TIBr 

sample 45-S10. 

Based on these initial results, which demonstrated the feasibility of bulk transport 
imaging in TIBr, further studies were performed to investigate the effects of doping and 
spatial variations. 
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V. THALLIUM BROMIDE SAMPLES ZRC0937-07 AND 

ZRC0937-08 


A. CRYSTAL GROWTH 

RMD provided slices from two different TIBr crystals doped with different 
impurities. RMD grew these crystals to investigate the efficiency of zone refining in 
removing various impurities. A synthesized 5N purity TIBr compound was purchased 
from a commercial vendor and used as the starting material. Crystal ZRC0937-07 was 
doped 1% atomic with 5N purity beads of NaBr, AlBr 3 , and AgBr. This crystal will be 
referred to as “Na/Al/Ag doped.” Crystal ZRC0937-08 was doped 1% atomic with 5N 
purity beads of CuBr 2 , FeBr 2 , and ZnBr 2 . This crystal will be referred to as “Cu/Fe/Zn 
doped.” Each crystal was sealed in a 10-inm diameter ampule under HBr atmosphere. 
The crystals were zone refined in a horizontal configuration for 100 passes using a 
travelling molten zone method. Crystal growth was done as a last step at a slow speed to 
promote single crystal growth. The grown crystals were cut into ~1 -mm thick slices using 
a diamond wire saw. The slices were grinded with a 30 microns alumina grit slurry to 
remove surface damages due to the wire saw cutting. They were then polished with a 
0.9 microns alumina grit slurry. The slices were etched with 2% Br in methanol. The 
grown crystals in their sealed ampules are shown in Figure 23. 
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(a) Na/AI/Ag doped TIBr crystal (ZRC0937-07) 



E 


Tail 


Head 


(b) Cu/Fe/Zn doped TIBr crystal (ZRC0937-08) 

A B C D E 



Figure 23. (a) Na/AI/Ag doped TIBr crystal (ZRC0937-07) and 

(b) Cu/Fe/Zn doped TIBr crystal (ZRC0937-08). 

The head region of each crystal contains fewer impurities 
than the tail region. 

Due to the segregation of impurities during zone refining, the head region of each 
crystal contains lower impurity concentrations than the tail region. RMD provided slices 
from different regions of the crystal for analysis. For the Na/AI/Ag doped crystal, RMD 
provided slices A and D, as shown in Figure 23(a). For the Cu/Fe/Zn doped sample, 
RMD provided slices A, C, and D, as shown in Figure 23(b). 

B. GLOW DISCHARGE MASS SPECTROMETRY 

RMD had EAG perform GDMS on TIBr slices near to the ones provided here for 
analysis. Figures 24(a) and (b) show the GDMS results for the Na/AI/Ag doped TIBr 
crystal and the Cu/Fe/Zn doped TIBr crystal, respectively. 
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Figure 24. GDMS results for (a) Na/Al/Ag doped TIBr and 
(b) Cu/Fe/Zn doped TIBr crystal. 


Figure 24(a) shows that concentrations of the dopants Na, Al, and Ag decrease 
from slice D to A (tail to head) in the Na/Al/Ag doped TIBr crystal. The majority of the 
impurities in this crystal also decrease in concentration from slice D to A. Figure 24(b) 
shows that the concentrations of dopants Cu, Fe, and Zn in the Cu/Fe/Zn doped TIBr 
crystal also decrease from slice D to A. The majority of the impurities in this crystal, 
however, increases in concentration or remain the same from slices D to A. The GDMS 
data indicates that these crystals, although zone refined, still have high concentrations of 
both dopants and impurities. 
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c. 


CATHODOLUMINESCENCE 


Room temperature (295 K) CL mapping was perfonned on slices A and D of the 
Na/Al/Ag doped TIBr crystal and on slices A and C of the Cu/Fe/Zn doped TIBr samples, 
as shown in Figure 25. The SEM operating parameters were 20.0 kV, lxlO' 9 A, and 
2000x magnification. The CL system was operated in panchromatic mode to capture 
emitted light at all wavelengths. The samples with higher impurity concentrations (slice 
D of Na/Al/Ag doped TIBr and slice C of Cu/Fe/Zn doped TIBr) show intensity 
variations on a scale of ~10 pm. 






Figure 25. CL map at 295 K of (a) slice A and (b) slice D of Na/Al/Ag doped TIBr 
and (c) slice A and (d) slice C of Cu/Fe/Zn doped TIBr. 

The room temperature luminescence patterns differ for slices A and C of the 
Cu/Fe/Zn doped TIBr, as seen in Figure 25(c) and (d). Slice A exhibits nearly uniform 
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intensity, whereas slice C shows spatial variations in the intensity on a scale of ~ 10 pm. 
CL mapping was repeated on these two slices at 11 K with SEM operating parameters of 
20.0 kV, 3 X 10" 10 A, and 2000 x magnification. The CL system was set to 460 nm 
(2.70 eV) to capture band edge luminescence. The CL maps are shown in Figure 26. 
Both slices now show intensity variations on a scale of ~ 10 pm. 

Slice A of the Cu/Fe/Zn doped sample showed unifonn intensity at room 
temperature, as seen in Figure 25(c), but spatial variation at 11 K, as seen in Figure 26(a). 
It is likely that slice A is actually non-uniform in material composition, but that the 
intensity of light emitted at 295 K is simply too low to reflect this non-unifonnity. It is 
also possible that slice A has regions that are more uniform than others, and that the CL 
map at 295 K was taken from one of these more uniform regions. 



Figure 26. CL map for X = 460 nm at 11 K of (a) slices A and 
(b) slice C of Cu/Fe/Zn doped TIBr. 


CL spectra were measured for the Na/Al/Ag and Cu/Fe/Zn doped TIBr samples 
at 10 K or 11 K. The SEM operating parameters were 20.0 kV, 3 x 10" 10 A, and 
2000 x magnification. For the CL system, the slits were set to 1.5 mm, and the step size 
was 1 nm. Figure 27(a-d) shows the CL spectra. All four spectra show the direct band 
gap at ~ 410 nm (3.02 eV) and the indirect band gap at ~ 460 nm (2.70 eV). 
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D. PHOTOLUMINESCENCE 


PL was also performed on the Na/Al/Ag and Cu/Fe/Zn doped samples at 
11 K using the PL system at the University of California, Berkeley. Figure 27(e-f) shows 
the CL spectra. All four spectra show the direct band gap at ~ 415 nm (3.00 eV) and the 
indirect gap at ~ 475 nm (2.61 eV). One sees that the midgap defect luminescence is 
more prominent in the CL spectra. This is most likely due to the larger energy density 
created by electron beam excitation, resulting in more phonon-assisted transitions. This 
should be a topic of further study. 


Cathodoluminescence 



Photoluminescence 



300 350 400 450 500 550 600 650 700 750 800 850 
Wavelength (nm) 


Figure 27. (a-d) CL and (e-f) PL spectra for Na/Al/Ag doped TIBr and 

Cu/Fe/Zn doped TIBr. 
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E. 


TRANSPORT IMAGING 


1. Effects of Doping 

The effect of overall doping concentrations on transport properties at 5 K was 
studied in the Na/Al/Ag doped and Cu/Fe/Zn doped TIBr crystal. First, single locations 
on each slice from each crystal were compared. Figure 28(a) shows a comparison of the 
averaged intensity distributions for slices A and D of the Na/Al/Ag doped TIBr crystal. 
Figure 28(b) shows a similar comparison for slices A and C of the Cu/Fe/Zn doped TIBr 
crystal. In each case, the intensity distribution is broader for slice A, the slice with the 
lower dopant concentration. This broadening indicates a larger /ur product, most likely 
associated with longer carrier lifetimes in the samples with lower impurity 
concentrations. 




Figure 28. Comparison of the normalized mean recombination luminescence 

intensity distribution for (a) slices A and D of the Na/Al/Ag doped TIBr crystal, 
and (b) slices A and C of the Cu/Fe/Zn doped TIBr crystal. Note that the figures 
are plotted on a logarithmic scale to highlight small changes in the distributions. 

Figure 29 shows the data and best fit lines from least square fitting to the model 
for slices A and D of the Na/Al/Ag doped TIBr crystal and for slices A and C of the 
Cu/Fe/Zn doped TIBr crystal. The first data point used for each fit is 4 pm from the 
source in order to avoid effects due to the generation volume created by the SEM beam. 
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Figure 29. Data and least squares fitting to the model for (a) (a) slices A and D 
of the Na/Al/Ag doped TIBr crystal, and (b) slices A and C of 
the Cu/Fe/Zn doped TIBr crystal. 


Table 2 shows the averaged ambipolar diffusion lengths and ambipolar /uz product 
for the four slices shown in Figure 29. The /uz product is higher in the samples with the 
lower impurity levels. This result is consistent with reported improvements in the 
electron /uz product and spectroscopic performance with improved material purification, 
although our measurements for the /uz product are dominated by the hole mobility, given 
the ambipolar nature of the diffusion [39], [40], [58], [59]. 


Table 2. GDMS dopant concentration, diffusion length, and /uz product at 5 K. 


GDMS Dopant Concentration (cm' 3 ) 


TIBr Crystal 

Slice 

Na 

A1 

Diffusion Length, L d 
Ag 

(pm) 

//r Product 

(cm 2 /V) 


A 

6.5xl0 16 

8.4xl0 15 

5.1xl0 16 

21 

l.OxlO' 2 

Na/Al/Ag doped 








D 

4.0xl0 18 

4.2xl0 16 

2.6xl0 17 

12 

3.3xl0' 3 


A 

<7.1xl0 14a 

4.1xl0 15 

< 6.9xl0 14a 

18 

7.7xlO' J 

Cu/Fe/Zn doped 








C 

6.4xl0 15 

9.0xl0 15 

<6.9xl0 14a 

9.0 

1.9xl0' 3 


a Below GDMS detection limit for given element 
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The data in Table 2 give a snapshot of the transport properties at one location on 
each slice and provide a gross estimate of how the transport properties change at different 
locations down the crystal ingot. The next step is to probe multiple locations on each 
slice to investigate how the uniformity of the /uz product varies within a slice. 

Multiple random locations on each slice were sampled to investigate the 
uniformity of the juz product across each slice. Figure 30(a) shows the intensity 
distributions for the Na/Al/Ag doped TIBr sample for slices A and D at random locations 
LI and L2. Figure 30(b) shows the intensity distributions for the Cu/Fe/Zn doped TIBr 
sample for slices A, C, and D at random locations LI, L2, and L3. 
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Figure 30. Intensity distributions on different slices (A, C, or D) for random 

locations (LI, L2, or L3) for (a) Na/Al/Ag doped and (b) Cu/Fe/Zn doped 

TIBr crystals at 5 K. 
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Figure 31 shows the extracted diffusion length for each slice from least squares 
fitting of the data to the model. 
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Figure 31. Diffusion lengths from different slices (A, C, or D) for random locations 
(LI, L2, and/or L3) for (a) Na/Al/Ag doped and (b) Cu/Fe/Zn doped 

TIBr crystals at 5 K. 


The data from Figure 31 are re-plotted in Figure 32 to show the mean and one 
standard deviation of the diffusion length. 
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Figure 32. Diffusion lengths and one standard deviation from different slices 
(A, C, or D) for random locations (LI, L2, or L3) for 
(a) Na/Al/Ag doped and (b) Cu/Fe/Zn doped TIBr crystals at 5 K. 
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The data shown in Figures 31 and 32 are summarized in Table 3. The associated 
/UT products and one standard deviation are also shown in the table. 


Table 3. Diffusion length and //r product at different locations on different slices of doped 

TIBr crystals at 5 K. 


TIBr Crystal 

Slice 

Location 

Diffusion Length, L d 

(pm) 

L d ± 1(7 

(pm) 

//rProduct 

(cm 2 /V) 

//r+ la 

(cm 2 /V) 



LI 

27 


1.7 X 10 2 



A 



26+2 


(1.5+0.3)x 10“ 2 

Na/Al/Ag doped 


L2 

24 


1.3xl0' 2 










LI 

7.3 


1.2x1 O' 3 



D 



9.6+3.3 


(2.3+1.4)xl0 -3 



L2 

12 


3.3x 10“ 3 




LI 

8.5 


1.7 X 10“ 3 



A 

L2 

12 

9.6+2.5 

3.6 X 10“ 3 

(2.2+1,2)x 10“ 3 



L3 

7.8 


1.4xl0 -2 




LI 

9.0 


1.9 X 10“ 3 


Cu/Fe/Zn doped 








C 

L2 

7.0 

11+6 

l.lxlO' 3 

(3.4+3.4)x 10' 3 



L3 

18 


7.3xl0' 3 




LI 

5.5 


7.1xl0' 4 



D 



5.4+0.2 


(6.8+0.4)xl0 -4 



L2 

5.3 


6.5xl0' 2 



We can draw several conclusions from the data in Table 3. For the Na/Al/Ag 
doped crystal, as we move down the crystal ingot from slice D to A (tail to head), the 
transport properties generally improve, consistent with overall reduction in impurity 
concentration due to zone refining. Within a slice, however, a comparison of the random 
locations shows that there is still significant variation in the transport properties across 
each slice. For the Cu/Fe/Zn doped sample, as we move down the crystal ingot from 
slice D to A, the transport properties do not strictly follow the expected trend of 
improvement. Instead, on average, slice C showed better transport properties than slice 
A. We can understand this variation by examining the transport properties of the random 
locations within slice C. Location 3 (L3) on slice C had a diffusion length double that of 
the other two locations on the slice. Most likely, L3 on slice C was a localized region of 
low dopant or impurity concentration. This again indicates significant variation in the 
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transport properties, most likely due to uneven distribution of dopants and impurities, 
both down the crystal ingot and across each slice. 

The CL maps of TIBr in Figures 25 and 26 showed spatial variations in the 
intensity, indicating an uneven distribution of impurities or defects. The TOF-SIMS 
imaging in Figures 17 and 18 also depicted variations in the impurity distribution. 
Transport imaging on multiple locations on each slice as shown in Figure 31 and 32 and 
summarized in Table 3 indicates that variations in these impurities or defects affect 
transport properties in the semiconductor. The next step is to use transport imaging to 
probe the scale of these variations. 

2. Spatial Variation in Diffusion Length and //r Product 

The ambipolar diffusion lengths and ambipolar //r products reported in Tables 2 
and 3 are average values over 40 pin regions of the TIBr slices at 5 K. By taking 
advantage of the resolution of the imaging system, however, we can obtain better spatial 
resolution. Transport imaging provides the ability to probe spatial variations in the 
diffusion length and //r products at a resolution of 2 pm [60]. This micro-scale technique 
is demonstrated on slice C of the Cu/Fe/Zn doped crystal. 

Figure 33(a) shows an optical image of a 60 pm x 85 pm SEM area luminescence 
scan at a random location on slice C of the Cu/Fe/Zn doped crystal at 5 K. The SEM 
operating parameters were 20.0 keV, 3 X 10‘ 10 A, 2000x magnification. This image shows 
significant spatial variation on a scale of ~10 pm. The image was acquired by exposing 
the sample to the SEM beam and recording the resulting luminescence on the CCD 
camera with an exposure time of 600 s. This image is similar to a CL map, but slightly 
different in that the luminescence information is obtained during area SEM excitation, 
rather than point by point. The location of a 40 pm line scan is indicated in Figure 33(a) 
by the red line. The line scan image is shown in Figure 33(b). This image was acquired 
with an exposure time of 90 s. The line scan is divided into twenty sections, each 2 pm 
in length. 
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(a) 



Optical luminescence 
from SEM picture mode 



Optical luminescence 
from SEM line scan 


Figure 33. (a) SEM area luminescence scan of 60 pm x 85 pm region of slice C of the 

Cu/Fe/Zn doped TIBr crystal. Intensity variations are visible on a scale of ~ 10 
pm. The location of a 40 pm line scan is indicated by the white line, (b) SEM 
line scan, divided into 2 pm segments. After [60]. 


The ambipolar diffusion length is determined from the mean intensity distribution 
over each 2-pm section. The results are shown in Figure 34. The best-fit ambipolar 
diffusion lengths range from 4.6 pm to 11.2 pm, with a mean of 7.1 pm and standard 
deviation of 2.0 pm. The corresponding pn products range from 4.9x 10‘ 4 cm 2 /V to 
2.9 X 10‘ 3 cm 2 /V, with a mean of 1.3><10‘ 3 cnr/V and standard deviation of 6.8 X 10‘ 4 
cm /V. These results demonstrate that, for these materials, CL luminescence variations 
are indicative of material variations that affect the ambipolar /jt product [60]. 
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Figure 34. (b) Diffusion lengths are averaged over 2 |xm sections of the line scan. 

The diffusion lengths vary between 4.6 jam and 11.2 jam, 
with a mean of 7.1 jam and standard deviation of 2.0 jam. After [60]. 

Crystal growers could use spatially resolved transport imaging as a micro-analysis 
technique to assess the quality of a crystal, without the need for contacts or device 
fabrication. Figure 35 shows a comparison of transport imaging for TIBr at 5 K and 72- 
type GaAs at 295 K. The figure shows that TIBr has a higher diffusion length, but that 
the crystal contains regions of nonuniformity that will likely affect detector performance. 
In contrast, although the GaAs sample has a lower diffusion length, the material is more 
uniform. Transport imaging can therefore provide rapid, yet still quantitative, assessment 
of how variations in crystal growth affect material uniformity. Further investigation may 
also provide a way to correlate the variations observed in transport imaging with detector 
resolution. 
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Figure 35. Variations observed in diffusion length in TIBr and GaAs using 

spatially resolved transport imaging. 
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VI. DEFECT STUDY IN SELENIUM AND LEAD-DOPED 

THALLIUM BROMIDE 


A. DEFECT LEVELS 

The level of numerous defect states in the bandgaps of Ge, Si, and GaAs are well 
known. Figure 36 shows a typical impurity level diagram for Ge. Diagrams like this 
have emerged as a result of over 40 years of research in Ge crystal growth, diffusion 
behavior, and materials characterization. For TIBr, however, theoretical and experimental 
research on the location of impurity levels is just now beginning to emerge [46], [61]. 
Cathodoluminescence in TIBr doped with Se and Pb can contribute to an understanding 
of the location of these impurity levels. 
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Figure 36. Impurity levels in Ge. The ionization energies are labeled in meV. 

Levels above the gap center (E,) are donors, unless denoted as an acceptor (A). 
Levels below the gap center are acceptors, unless denoted as a donor (D). 

After [47]. 


B. CRYSTAL GROWTH 

TIBr crystals doped with either Se or Pb were investigated in this study. The 
crystals were grown by RMD using the Vertical Bridgeman technique. One crystal was 
nominally doped with 100 ppm T^Se, and the other crystal was nominally doped with 
100 ppm PbBr 2 . The crystals are shown in Figure 37. The figure shows the grown 
crystals and the approximate locations of the slices, indicated by red lines. For each 
crystal, slice Cl was used for our analysis. 
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(a) Se doped TIBr crystal (VB-D01) 
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(a) Pb doped TIBr crystal (VB-D03) 
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Figure 37. TIBr crystals doped with either Se or Pb. Approximate locations of slices 
are shown in red. (a) For the Se doped TIBr crystal, slice C3 is shown, 
although slice Cl was used in this study, (b) For the Pb doped TIBr crystal, 

slice Cl was used. 


C. GLOW DISCHARGE MASS SPECTROMETRY 


GDMS was performed by EAG on the Se and Pb doped TIBr crystals. For the Se 
doped crystal, GDMS was perfonned on slice Cl, the same slices used in this study. For 
the Pb doped crystal, GDMS was performed on slice C2, the slice next to the one used in 
this study. The GDMS results are shown in Figure 38. The Se doped crystal shows high 
levels of both Se and Cl. The Pb doped crystal shows high levels of Pb, C, N, and Cl. 
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Figure 38. GDMS results for (a) Se doped TIBr and (b) Pb doped TIBr crystal. 
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D. 


CATHODOLUMINESCENCE 


Room temperature (295 K) and low temperature (5 K) CL spectroscopy were 
performed on the Se and Pb doped TIBr slices to confirm the presence of a CL signature. 
The SEM operating parameters were 1><10‘ 9 A and 500 x magnification. The spectra for 
the Se and Pb doped TIBr slices are shown in Figure 39. The spectrum for the Pb doped 
sample in Figure 39(b) includes Gaussian fits to peak C and D. 




Figure 39. CL spectra at 295 K and 5 K for (a) Se doped TIBr and (b) Pb doped TIBr. A 

Gaussian fit to peak D is shown as a dashed line. 


Figure 40 compares the spectra for the Se and Pb doped samples at 5 K. The Se 
doped sample shows peaks at 3.01 eV and 2.67 eV (A and B) due to direct and indirect 
excitons [53]. The spectrum also shows a defect peak at 2.43 eV (C) and a wide defect 
band at 1.87 eV (D). The Pb doped sample also shows the direct and indirect peaks at 
3.02 eV and 2.69 eV (A and B). The peak at 2.43 eV observed in the Se doped sample is 
absent. The Pb doped sample shows additional peaks at 2.23 eV (C) and 2.12 eV that 
were not observed in the Se doped sample. 
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Figure 40. CL spectra at 5 K for (a) Se doped TIBr and (b) Pb doped TIBr. 


E. DEEP LEVELS IN THALLIUM BROMIDE 

The CL spectra for Se and Pb doped TIBr both show mid-gap defect peaks. The 
CL spectroscopy can be used in conjunction with other spectroscopic techniques to 
identify the dopant-related deep level impurity traps that cause these peaks and to begin 
to piece together a defect level diagram for TIBr similar to the one shown in Figure 36 for 
Ge. Theoretical work using Density Functional Theory (DFT) has been done to identify 
deep level impurity sites in the bandgap of TIBr [46]. To the author’s knowledge, 
however, there are currently no published band diagrams based on experimental work 
that show the location of dopant-related deep levels in TIBr. 

One technique for studying deep levels in semi-insulating semiconductors is 
photo-induced conductivity transient spectroscopy (PICTS). PICTS provides information 
on the trap activation energy, the trap capture cross section, and the trap type (electron or 
hole). Together, PICTS and CL can be used to identify the location of deep levels in a 
semiconductor [62]. Our collaborators at the University of California, Berkeley, have 
performed PICTS on our samples of Se and Pb doped TIBr. By combining PICTS results 
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with the CL spectroscopy, the first band picture of dopant-related deep levels for Se and 
Pb has been created [61]. The results are briefly summarized here. 

In the PICTS experimental set up at Berkeley, a high-brightness LED peaked at 
450 nrn (2.76 eV) is used as a photoexcitation source. A bias between 15-50 V is applied 
to the sample, and the light is pulsed on for 350 ms and off for 150 ms. The transient 
photoresponse is measured as a function of temperature over the range of 80-250 K in 
0.5 K increments. PICTS spectra are generated, as shown in Figure 41(a). Peaks in the 
spectra denote detected traps. The traps are plotted on an Arrhenius plot, as shown in 
Figure 41(b). The slope of the line on the Arrhenius plots gives the trap activation 
energy, and the y-intercept provides the trap capture cross section. For Se doped TIBr, 
two traps (A and B) were identified. For Pb doped TIBr, three traps (C, D, and E) were 
identified. 



Figure 41. (a) PICTS spectrum and (b) Arrhenius plot. After [61]. 


For Se and Pb doped TIBr, the measured activation energy, capture cross section, 
and trap type are summarized in Table 4. 


Table 4. Activation energies, cross section, and trap type for identified deep levels. 

After [61] 


Peak 

Sample Dopant 

E activation 

feV] 

Capture Cross Section 
|em 2 l 

Tentative Trap Type 

A 

Se 

0.35 

1.9 x 10" 12 

electron 

B 

Se 

0.45 

6.6 x 10 12 

hole 

C 

Pb 

0.11 

1.1 x 10' 16 

electron 

D 

Pb 

0.45 

3.9* 

hole 

E 

Pb 

0.75 

2.4 x 10 3 * 

hole 


* Denotes unphysically large capture cross sections. 

65 













From the CL spectra, an energy offset, AE ind , is measured between the defect peak 
and the indirect band edge at 2.68 eV. Transitions between the levels identified in the 
PICTS analysis are then matched to AE in d to create the band picture. Table 5 shows the 
dopant, CL peak, AE inc /, and matching transitions from PICTS. For the Se doped TIBr 
sample, the CL peak at 2.43 eV does not match a transition identified with PICTS. The 
CL peak at 1.86 eV, however, matches an intradefect transition between the electron trap 
at 0.35 eV below the conduction band edge (denoted E c — 0.35 eV) and the hole trap at 
0.45 eV above the valence band edge (denoted E v + 0.45 eV). For the Pb doped TIBr 
sample, the CL peak at 2.23 eV matches a transition from the conduction band edge to 
the hole trap at E v + 0.45 eV. The CL peak at 2.12 eV matches an intradefect transition 
from the electron trap at E c — 0.11 eV to the hole trap at E v + 0.45 eV. The PICTS level at 
E v + 0.75 does not match a CL peak. 


Table 5. Comparison of CL and PICTS data. After [61] 


Dopant 

CL Peak [eV] 

AEind 

PICTS [eV] 

Se 

2.43 

0.25 

- 

Se 

1.86 

0.82 

0.35 (A) + 0.45 (B) 

Pb 

2.23 

0.45 

0.45 (D) 

Pb 

2.12 

0.56 

0.45 (D) + 0.11(C) 

Pb 

- 

- 

0.75 (E) 


Together, the CL and PICTS data result in the first empirical energy diagram 
showing deep level impurity states related to Pb and Se dopants in TIBr. The diagram is 
shown in Figure 42. 
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Figure 42. Energy diagram showing mid-gap defect states related to 
Pb and Se dopants in TIBr. Diagram not to scale. After [61] 


F. TRANSPORT IMAGING 


Transport imaging was performed on the Se and Pb doped TIBr samples to 
determine the effect of doping on the ambipolar diffusion length. Three random locations 
on each sample were chosen. Figure 43 shows the intensity distributions for each 
location on each sample at 295 K and 5 K. For each sample, at each temperature, there is 
little variation in the intensity distributions across the different locations. This indicates 
good material uniformity. Also, the distributions broaden at low temperatures, indicating 
increased carrier mobilities and lifetimes at lower temperatures. 




Intensity distributions for three random locations at 295 K and 
5 K for (a) Se doped TIBr and (b) Pb doped TIBr. 
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Figure 43. 



















A comparison of the intensity distributions for each sample at 5 K is shown in 
Figure 44. The intensity distributions are significantly broader for the Se doped sample. 
This broadening is consistent with the GDMS data shown in Figure 38, which indicates 
that the Se doped TIBr contains fewer impurities than the Pb doped sample. However, 
since each sample still has high impurity levels, Pb alone may not be entirely responsible 
for the decreased intensity distributions observed in Figure 44. Lead may have the 
dominant effect, but other impurities are also likely to affect the transport properties. 
Further investigation with samples with a range of dopant types and concentrations would 
help to determine the cause of the diffusion length decrease in the Pb doped TIBr sample. 



Figure 44. Intensity distributions for three random locations at 5 K for 

Se doped TIBr and Pb doped TIBr. 


The mathematical fits differ for the Se and Pb doped samples. The fits are shown 
for each sample in Figure 45. For each sample, the lower panel is shown on a 
logarithmic scale and highlights the differences between the fits. The fit for the Se doped 
TIBr sample follows the data for the full data range, out to 50 pm. The fit for the Pb 
doped TIBr sample, however, shows systematic deviation from the data beginning at 
approximately 15 pm. Further studies using a range of Pb doped samples may help to 
understand this deviation. 
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Figure 45. Least squares fitting to the data at 5 K for (a) Se doped TIBr 

and (b) Pb doped TIBr. 

The diffusion lengths from least squares fitting to the model and the calculated /lit 
products for each sample at each location at 5 K are shown in Table 6. 


Table 6. Diffusion length and //rproduct at random locations 

on doped TIBr crystals at 5 K. 


TIBr Crystal 

Location 

Diffusion Length, L d 

(pm) 

L d ± 1 a 

(pm) 

fir Product 

(cm 2 /V) 

fir± 1 a 

(cm 2 /V) 


1 

14 


4.6 X 10" 3 


Se doped 

2 

12 

13±1 

3.2><10" 3 

(3.8+0.7)xl0 -3 


3 

13 


3.7 X 10" 3 



1 

6.7 


l.OxlO -3 


Pb doped 

2 

6.2 

6.5+0.3 

8.9xl0' 4 

(8.9+0.6)xl0 -4 


3 

6.7 


l.OxlO -3 



Figure 46 shows the spatial variation in the diffusion lengths along 16 pm 
sections of the Se and Pb doped TIBr, averaged over 2 pm sections of the line scan 
images at 5 K. The materials are the most uniform of any of the TIBr materials 
investigated to date in terms of their diffusion length and /tr product. Transport imaging 
is again able to provide a rapid assessment of both the magnitude of the transport 
properties and the material uniformity in these bulk samples 
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Figure 46. Spatial variation in the diffusion lengths along 16 pm sections of the 
Se and Pb doped TIBr, averaged over 2 pm sections of the 
line scan images, at 5 K. 
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VII. INITIAL INVESTIGATION IN BULK CADMIUM ZINC 

TELLURIDE 

A. OVERVIEW OF CADMIUM ZINC TELLURIDE 

Cadmium zinc telluride (CZT) is one of the most promising compound 
semiconductors under development as a room temperature replacement for HPGe for 
nuclear radiation detection [24], [56], [63] - [65], Like TIBr, CZT has a high atomic 
number (Zed = 48, Z Zn = 30, Z Te = 52) and high density (5.78 g/cm ), resulting in high x- 
ray and gamma photon attenuation. The wide bandgap of CZT (1.57 eV) leads to high 
resistivity (~10 10 Q-cm), which allows for room temperature operation without excessive 
thermal generation of carriers and dark current [66], [67]. Unlike TIBr, CZT exhibits no 
room temperature polarization effects, which allows for long-term operation without 
crystal degradation [24], [56]. In terms of spectroscopy, energy resolution < 1% FWHM 
has been achieved with CZT for Cs (662 keV) [64], Table 7 summarizes select 
properties of CZT. 


Table 7. 

Properties of CZT. After [21] 

Property 

Value 

Atomic Number 

Z ci \ 48, Z Zn : 30, Z Te : 52 

Density 

5.78 g/cm 3 

Bandgap 

1.57 eV 

Resistivity 

3><10 10 Q-cm 

Electron Mobility 

1000 cnf/V-s 

Hole Mobility 

120 cm 2 /V-s 

Electron /it Product 

4xlO _3 cm 2 /V 

Hole /rr Product 

1.2xl0" 6 cm 2 /V 

Melting Point 

1092-1295 °C 

Crystal Structure 

Cubic (Zinc Blende) 
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Despite its many positive attributes, several obstacles have limited the widespread 
use of CZT [68]. Large volume, inexpensive, high quality crystals have been difficult to 
grow [24], CZT suffers from incomplete charge collection, particularly for holes, which 
results in spectral broadening [56], [69]. To overcome this deficiency, CZT detectors are 
often operated as “electron-only” devices that only collect electrons [69], Electron 
charge collection, however, is still hampered by the presence of Te inclusions. Te 
inclusions are non-stochiometric defects that form during the melt growth of the crystal, 
with a typical size of 1-2 pm. Their presence and non-uniform distribution in CZT 
degrades energy resolution [68]. 

B. CATHODOLUMINESCENCE 

CL mapping and spectroscopy was performed on CZT samples provided by 
Redlen Technologies in order to determine the CL signature. Two samples, CZT(STl) 
and CZT(372), were investigated. The difference between the two samples is that 
CZT(STl) has been annealed, but CZT(372) has not. 

CL spectroscopy was performed at 295 K and 5 K. Lor CL spectroscopy, the 
SEM operating parameters were 20.0 keV, 3 X 10' 10 A, and 2000 x magnification. The slits 
were set to 0.5 mm (295 K) or 0.1 mm (5 K), and the monochromator step size was 1 mn. 

The CL spectra for each sample are shown in figure 47. for the spectra at 295 K, 
the peak at 1.56 eV is the band edge of CZT. for the spectra at 5 K, the peak at 1.64 is 
the band edge of CZT. The peak at -1.61 eV is likely the second LO-phonon (21 meV) 
replica [70]. The peak at - 1.47 eV is likely defect luminescence from complex deep 
level defects involving Cd vacancies and impurities [71]. 
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Figure 47. CL spectra at 295 K and 5 K for (a) CZT(STl) and (b) CZT(372). 


A comparison of the CL spectra for each sample at 5 K is shown in Figure 48. 
The most noticeable feature is a slight decrease in the defect luminescence at -840 nm 
(1.48 eV) in the unannealed sample, CZT(372). 
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CL mapping was performed on each spectrum at 5K at 755 nm (1.64 eV) to 
capture the band edge luminescence. The SEM operating parameters were 20.0 keV, 
3 X 10" 10 A, and 1000x magnification, and the slits of the CL system were set to 1 mm. 
The CL maps are shown in Figure 49. The CL map of CZT(STl), as shown in Figure 
49(a), exhibits areas of decreased intensity in approximately triangular shaped regions, 
with sizes of 1-5 pm. The shape and size of these regions are consistent with the 
reported shape and size of Te inclusions in CZT [50], The CL map of CZT(372), as 
shown in Figure 49(b), exhibits black dots up to ~1 pm each along a diagonal line. These 
dots may represent Te inclusions that have migrated to a defect, such as a grain boundary, 
during cooling after the growth of a CZT ingot. Such Te inclusion along a defect are said 
to “decorate” the defect [48]. Further correlated studies with CL and EDX would be 
necessary to positively identify these areas of reduced intensity as Te. 



Figure 49. CL mapping at 5 K for (a) CZT(STl) and (b) CZT(372). 


C. TRANSPORT IMAGING 
1. Temperature Studies 

Transport imaging was performed on CZT(STl) and CZT(362) at 295 K to 
investigate their room transport properties. The SEM operating parameters were 10.0 
keV, 6 X 10‘ 10 A, and 500 x magnification. Figure 50 shows the intensity distributions for 
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each sample. For each sample, two random locations, labeled LI and L2, were sampled. 
Each location represents a single field of view from the SEM in picture mode, and these 
locations on the crystal could be separated from each other by a distance of 0.5 mm or 
more. At each location, two separate line scan images were recorded, at two closely 
spaced sub-locations. First, one line scan image was recorded. Then, without moving the 
sample or changing the field of view, the SEM beam was moved parallel to the first line 
scan image, and a second line scan image was recorded. Each line scan image is labeled 
with the pixel sub-location in the field of view. In Figure 50(b), for example, locations 
LI and L2 on sample CZT(372) represent two different random locations, separated by 
0.5 mm or more. The separation between sub-locations LI x624 and LI x250, however, 
is only 624 - 250 = 374 pixels. Based on the system resolution of 0.4 pm/pixel, these 
two sampled sub-locations are only separated by ~ 150 pm. The line scan images for 
these two sampled sub-locations are overlaid in the inset in Figure 50(b). Using this 
scheme, we are able to investigate possible gross and fine changes in transport properties 
across the crystal. As shown in Figures 50(a) and (b), the intensity distributions for each 
sample are approximately the same for each location, indicating good material 
uniformity. 




Figure 50. Transport imaging at two locations, LI and L2, for samples (a) CZT(STl) 

and (b) CZT(372). Locations LI and L2 represent widely separated locations on 
the crystals, up to 0.5 mm or more. LI x492 and LI x391, for example, represent 
sub-locations that may only be separated by 150 pm or less. The inset shows an 
overlay of the two line scan images at sub-locations LI x250 and LI x624. 
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A comparison of the intensity distributions of the two samples at two random 
locations, LI and L2, at 295 K is shown in Figure 51. For CZT(STl), the distributions 
for the two sub-locations of LI (x492 and x391) are nearly identical. For CZT(372), the 
distributions for the two sub-locations of L2 (x275 and xl 13) are nearly identical as well. 
At first glance, when comparing the distributions for CZT(STl) to the distributions of 
CZT(372), it may appear that the distributions for CZT(STl) are more broad. Over the 
region of approximately 4-12 pm, however, the slopes of the lines are nearly identical. 
This indicates that the diffusion lengths are approximately equal in these materials. 



Figure 51. Comparison of intensity distribution at two different locations 

at 295 K for CZT(STl) and CZT(372). 

Transport imaging was repeated on samples CZT(STl) and CZT(372) at 5 K in 
order to detect any low temperature differences in diffusion lengths between the two 
materials. Figure 52 shows the intensity distributions at 5 K at two closely spaced sub¬ 
locations on CZT(STl) and two different locations on CZT(372). Similar to the results at 
295 K, the slopes of the distributions are nearly identical beginning ~ 4 pm from the 
source, after an initial difference in the shoulder regions of the curves. 
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Figure 52. Comparison of intensity distribution at 5 K for CZT(STl) and CZT(372). 

For least squares fitting to the model for CZT(STl) and CZT(372) at 295 K and 5 
K, data were used beginning at 4 pm from the source. The last data point used was 
varied between 16 pm, 22 pm, and 50 pm. Table 8 shows the calculated diffusion 
lengths at each temperature, for each sample, for each location, for each stopping point. 
The reported “endpoint average” diffusion length is an average value over all stopping 
points for each location. The reported “sample average” is an average value over all 
stopping points over all locations for each sample. The reported /uz product is based on 
this sample average. 
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Table 8. Ambipolar diffusion length and //rproduct for CZT(STl) and CZT(372) at 295 K and 5 K. 







Diffusion Length, L,i (pm) 








Data Used for Fitting 












Endpoint 

Sample 

/jt Product 

Temp. 

Sample 

Location 

4-12 pm 

4-15 pm 

4-22 pm 

4-50 pm 

Average, 

Average, 




Ld,End 

Ld,Sample i 

(cm'/V) 







LI x492 

2.67 

2.69 

2.69 

2.69 

2.69 




CZT(STl) 







2.68 + 0.01 

2.83x10^ 

295 K 


LI x391 

2.66 

2.67 

2.67 

2.67 

2.67 




CZT(372) 

L2 x275 

2.78 

2.77 

2.76 

2.75 

2.77 

2.72 + 0.05 

2.91xl0‘ 6 


L2 xl 13 

2.70 

2.68 

2.66 

2.66 

2.68 





LI x733 

3.08 

3.11 

3.14 

3.14 

3.12 




CZT(STl) 







3.05 + 0.08 

2.16x10“* 

5 K 


LI x810 

2.94 

2.97 

2.99 

2.99 

2.97 




CZT(372) 

LI x623 

3.14 

3.29 

3.36 

3.37 

3.29 

3.02 + 0.30 

2.1 lx 10~ 4 


L2 x702 

2.68 

2.74 

2.77 

2.77 

2.74 
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There are three conclusions to draw from the table. First, the stopping point used 
has little impact on the fit for the diffusion length. At 295 K for CZT(STl) LI x492, for 
example, the diffusion length is 2.67 pm based on data from 4-12 pm, and 2.69 pm 
based on data from 4-15 pm, 4-22 pm, and 4-50 pm. For all fitting in this dissertation, 
50 pm has therefore been used as the stopping point. Second, the two samples do not 
exhibit significant differences in diffusion length at a specific temperature, such as 295 
K. At 295 K for CST(STl) and CZT(372), for example, the sample average diffusion 
lengths are 2.68+0.01 pm and 2.72+0.05 pm, with only a slight difference between them. 
Third, the samples may exhibit detectable differences between temperatures. At 295 K 
for CZT(372) L2 x275, for example, the endpoint average diffusion length is 2.77 pm. 
At 5 K for CZT(372) Li x623, however, the endpoint average diffusion length is 
3.29 pm. Therefore, it may be possible to track changes of diffusion length with 
temperature. 

Table 8 showed that the end point used for data fitting has little impact on the 
calculated diffusion length. The start point, however, has a significant impact on the fit. 
Figure 53 shows fits for CZT(STl) LI x733 at 5 K for the data range 2-12 pm and 
4-12 pm. For 2-12 pm, the diffusion length is 4.71 pm, but for 4-12 pm, the diffusion 
length is 3.08 pm. As shown in Figure 53(a), the model cannot properly fit the shoulder 
region of the curve exhibited in this sample below ~ 4 pm. Further investigation is 
required to better understand this region right near the source. 
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Figure 53. Mathematical fits for diffusion length at 5 K for CZT(STl) LI x733 

for data from (a) 2-12 |um and (b) 4-12 jam. 


Table 8 showed that it may be possible to detect changes in diffusion length as a 
function of temperature. A temperature study was conducted on CZT(372) to investigate 
how diffusion length and /jt product varied over the range of 5 K to 60 K. Figure 54 
shows the normalized mean distributions for the temperature range. The SEM operating 
parameters were 10.0 keV, 6 X 10" 10 A, and 1000* magnification. 
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Figure 54. Normalized mean intensity distributions for CZT(372) for temperatures 

between 5 K and 60 K. 
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Figure 55 shows plots of temperature vs. diffusion length and temperature vs. /it 
product obtained by fitting the data in Figure 54 to the model. The /lit product decreases 
monotonically with increasing temperature over the range of 5-60 K. 



Figure 55. Temperature dependence of (a) diffusion length and (b) px product 

for CZT(372). 


2. Effects of Inclusions 

Detector-grade CZT is typically grown in Te-rich melts, which results in ingots 
with high concentrations of Te precipitates and Te inclusions. Te precipitates form 
during the cooling process, and their typical diameter is 10-50 mn. In contrast, Te 
inclusions are non-stochiometric defects that form during the crystal melt growth as a 
result of the capture of Te-rich melt droplets at the growth interface. Their typical size is 
1-2 pm, with sizes up to 100 pm reported, and they have a characteristic triangular 
shape. A dense dislocation field typically surrounds a Te inclusion. Te inclusions can 
also migrate during the cooling of the CZT ingot. When they encounter a structural 
defect, such as a grain boundary, their movement through the crystal stops, and the Te 
inclusions decorate the defect. Te inclusions and their associated dislocation field can 
degrade charge transport by trapping charge and distorting the internal electric field of 
the crystal. Non-uniform distribution of Te inclusions can result in spatial variation of 
charge transport and degraded energy resolution [50], [68]. Not all inclusions, however, 
affect charge transport. Te inclusions with diameters < 3 pm have been shown to have no 

effect on charge transport, up to concentrations of 3><10 cm' [68]. 
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Amman (2002) has proposed a model to explain the non-uniformity in electron 
transport observed in CZT crystals [26]. His model is shown schematically in Figure 56. 
In this model, one alpha particle (ai) passes through a region free of Te inclusions and 
creates charges that are trapped uniformly through the crystal. A second alpha particle 
( 012 ) passes through a region with a high concentration of inclusion and experiences non- 
uniform charge trapping. The two events will register different pulse heights in the 
detector, resulting in spectral broadening and degraded detector resolution. 


Alpha- 

particle 

beam 

diameter 


inclusion Region of increased 
v electron trapping 

A 



Alpha-particle char & e tapping 

interaction 
volume 


Figure 56. Proposed explanation for non-uniform charge trapping 

due to Te inclusions. After [26]. 


An initial transport imaging experiment has been performed to illustrate a 
potential approach to investigating the effects of Te inclusion on charge transport. Figure 
57(a) shows an optical image of CZT(STl) from the SEM in picture excitation mode at 
5 K. Four dark spots are visible, each ~ 1 pm in width, and two of them, labeled x733 
and x810, exhibit clear triangular shapes. Based on their sizes and shapes, these dark 
triangular spots are presumably Te inclusions. Figure 57(b) shows the normalized 
luminescence intensity from the same region of the crystal. The luminescence from the 
Te inclusions is very low. 
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Figure 57. (a) Luminescence from CZT(STl) from SEM operated in picture 

excitation mode, (b) Luminescence intensity from same region 
of CZT(STl) crystal. 

Figure 58(a) shows a slightly extended area of the same region of the CZT(STl) 
crystal. The red line indicates where the SEM line mode excitation will cross right over 
the Te inclusion labeled x733. An optical image of the line scan is shown in Figure 
58(b). The line scan data are analyzed over two regions. Region B-B' includes 160 pm 
of data. Region A-A' primarily includes the Te inclusion. The intensity distribution for 
each region is shown in Figure 58(c). Surprisingly, the intensity distribution is broader 
over the inclusion compared to the average distribution over the entire line scan. From 
the least squares fitting of the data from the two regions, the diffusion length is 3.07 pm 
over the entire line scan (B-B') but increases to 4.05 pm in the vicinity of the Te 
inclusion (A-A'). The process was repeated with the Te inclusion labeled x810, with 
similar results. The diffusion length was 2.99 pm over the entire line scan, and 4.15 pm 
over just the inclusion. 
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Figure 58. (a) Optical image from SEM picture excitation mode, showing location 

of line scan in red. (b) Optical image from SEM line excitation mode. The data 
are extracted over the majority of the line scan (B-B') and over just the Te 
inclusion (A-A'). (c) Luminescence intensity distributions over the two regions. 

Further investigation is necessary to gather more statistics. These results could be 
within the experimental error, or it could be that the Te inclusion is a sink for local point 
defects and the diffusion length is in fact locally increased in its vicinity. 

Using transport imaging, it may be possible to quantify the region of increased 
electron trapping associated with a Te inclusion. A series of line scan images could be 
taken, beginning at one side of the inclusion and continuing to the other side, as shown 
schematically in Figure 59. Each line scan could be analyzed twice, once over the entire 
line scan (B-B'), and a second time just over the region in line with the inclusion (A-A'). 
Far to the left of the inclusion, outside the region of increased electron transport, we 
would expect the two regions from each line scan to have approximately the same 
diffusion length. As we approach the region of influence, however, we would expect the 
diffusion behavior of the two line scan regions to diverge. On the right side of the 
inclusion, we would similarly expect to approach a point at which the diffusion behavior 
from the two regions converges again. One could further apply this technique 
systematically to inclusions of varying sizes to investigate the correlation between 
inclusion size and its influence on the region of increased electron trapping. 
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Figure 59. Schematic for investigating region of increased electron trapping 

associated with Te inclusions. 
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VIII. FUTURE WORK 


A. OVERVIEW OF FUTURE WORK 

In this work, we have applied transport imaging to bulk semiconductors to 
investigate spatial uniformity in the /lit product. Fitting the data to the mathematical 
model remains challenging in two regions. Close to the source, at distances less than 
~ 4 pm, a shoulder region often appears in the data, which the model cannot 
accommodate. At distances greater than ~15 pm, when intensity values reach < 1% of 
the maximum, systematic deviation often occurs between the data and the mathematical 
model. This is most likely due to statistical fluctuations associated with background 
levels in the CCD camera. Future work using transport imaging in bulk samples should 
focus on modifying the model to incorporate the physics at play in these two regions. 
Three areas of potential research are photon recycling, the generation volume, and a high- 
injection model. 

B. PHOTON RECYCLING 

Photon recycling, also known as self-absorption, occurs when a high-energy 
photon created by the recombination of an EHP is re-absorbed to create another EHP 
nearby. This new EHP can then recombine and re-emit a photon of lower energy. The 
net result is the production of lower energy photons of longer wavelength. Photon 
recycling has been studied extensively for a variety of photovoltaic and optoelectronic 
devices [72] - [79]. In double heterostructures of GaAs-GaAlAs, for example, photon 
recycling has been shown to increase the net radiative lifetime of carriers [79]. Of special 
note, Badescu (1993 and 1997) has developed a model to show how photon recycling 
influences the spatial distribution of excess minority carriers in n-type or />typc 
semiconductors under low injection conditions and results in increased efficiency in 
GaAs solar cells [74], [75]. 
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The effects of photon recycling in bulk semiconductors such as TIBr have not 
been investigated. Here, we report on preliminary results of photon recycling in CZT and 
TIBr using transport imaging with optical filters. 

In transport imaging, charge transport infonnation is obtained from the 
luminescence intensity distribution. To accurately model the diffusion behavior of the 
charge carriers in semi-insulating materials, this luminescence should only come from 
EHP recombination, not secondary emission due to photon recycling. The potential 
effect of photon recycling on the intensity distribution is illustrated schematically in 
Figure 60. In Figure 60(a), no photon recycling is present. EHPs created in the 
generation volume diffuse, and then recombine, emitting photons of various wavelengths, 
depending on the luminescence spectrum of the individual material. The intensity 
distribution reflects the luminescence due to the diffusion of EHPs alone. In Figure 
60(b), photon recycling occurs. EHPs created in the generation volume diffuse, and then 
recombine, emitting photons of different wavelengths. Some of the high energy photons, 
however, may create secondary photons of lower energy at points removed from the site 
of the EHP recombination, especially if midgap defects are present in the material. The 
overall effect is to broaden the intensity distribution. Now, the intensity distribution 
represents both diffusion and photon recycling. 




Figure 60. Effects of photon recycling in a semiconductor. 
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One way to investigate the role of photon recycling is to filter the light before it 
reaches the CCD array camera. The CL spectra for both CZT(372) and TlBr(Se) were 
observed to have high energy peaks followed by lower energy broad defect 
luminescence. Transport imaging was performed in these samples with short and long 
pass filters to separate the luminescence observed. 

If photon recycling is occurring in these samples, we would expect the following 
results. With the short pass filter, we should primarily detect the shorter wavelength, 
higher energy photons created by the recombination of EHPs. The intensity distribution 
should be narrower than the unfiltered distribution since we are not detecting secondary, 
lower energy photons created by photon recycling at points far removed from the 
generation volume. With the long pass filter, we should primarily detect the longer 
wavelength, lower energy photons, created both by the recombination of EHPs and by 
photon recycling at points far removed from the generation volume. The intensity 
distribution should actually be broader than the unfiltered spectrum, however, since we 
are lowering the overall recorded intensity but still capturing light at the same distance as 
the unfiltered distribution. 

Figures 61 shows the CL spectrum at 5 K for CZT(372) and the transmission 
curves for the optical filters used. The CZT(372) CL spectra at 5 K shows a direct gap 
peak at 756 nm (1.64 eV) and a broad defect peak at 839 nm (1.48 eV). As shown in 
Figures 61(a) and (b), a short pass filter at 800 nm was used and a long pass filter at 
780 nm was used. 



Wavelenth (nm) Wavelenth (nm) 


Figure 61. CL spectra at 5 K for CZT(372) and filter transmission spectra for (a) Short pass 

filters and (b) Long pass filters. 
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Figure 62 shows the intensity distribution from transport imaging of CZT(372). 
Figure 62(a) shows the intensity plotted in terms of counts recorded on the CCD array 
camera, and Figure 62(b) shows the same data plotted with normalized intensity. 



Position (fiin) 



Figure 62. Intensity distributions from transport imaging of CZT(372) at 5 K for unfiltered, 
short-pass filtered, and long-pass filtered conditions. The same data are shown 
with the y-axis plotted in (a) Intensity, as measured by counts recorded on the 
CCD array camera, and (b) Normalized scale, with each distribution normalized 

to its own maximum. 


The intensity distributions shown in Figure 62 for CZT(372) show the trend we 
would expect if photon recycling is occurring. The distribution with the short pass filter 
is narrower than the unfiltered distribution due to the removal of secondary, lower energy 
photons created by photon recycling at points far removed from the generation volume. 
The distribution with the long pass filter is broader than the unfiltered distribution due to 
removal of higher energy photons that reduce with intensity while still detecting lower 
energy photons created through photon recycling at points far removed from the 
generation volume. Although further investigation is necessary, this result suggests that 
future transport imaging work in CZT should be conducted with a short pass filter to 
remove the effects of photon recycling. 

Figure 63 shows the CL spectrum at 5 K for TlBr(Se) and the transmission curves 
for the optical filters used. The TlBr(Se) CL spectra at 5 K shows direct and indirect gap 
peaks at 412 nm (3.01 eV) and 465 nm (2.67 eV) and a broad defect peak at 664 mn 
(1.87 eV). As shown in Figure 60(a), two short pass filters at 550 nm and 600 nm were 
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used. The filters were used in tandem due to the increase in their transmission at longer 
wavelengths. Since the transmission increases at longer wavelengths, it was not possible 
to completely filter these longer wavelengths. As shown in Figure 60(b), a long pass 
filter at 515 mn was used. 
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Figure 63. CL spectra at 5 K for TlBr(Se) and filter transmission spectra for (a) Short pass 

filters and (b) Long pass filters 


Figure 64 shows the intensity distribution from transport imaging of TIBr. Out to 
~4 pm, the distributions follow the trend we would expect if photon recycling were 
present. This is not the case, however, beyond ~4 pm, where the distribution is broader 
with the short pass filters. This could be due to incomplete filtering of the longer 
wavelength photons with the short pass filter, allowing some photon recycling to occur 
and broadening the distribution. It could also suggest that carrier diffusion is the 
dominant process and that the short pass distribution simply reflects variations in the 
recombination ratios. Further investigation with different filters is necessary to 
understand these results. These initial results do not suggest that traditional photon 
recycling is dominant in TIBr. 
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Figure 64. Intensity distributions from transport imaging of TlBr(Se) at 5 K for unfiltered, 

short-pass filtered, and long-pass filtered conditions. 


The 2D diffusion model for charge transport in bulk semiconductors does not 
include a term to account for photon recycling. Many methods have been investigated 
for including this term [77] - [79]. One approach [78] is to include a radiative lifetime 
multiplicative factor, (j), whereby the minority carrier lifetime becomes 


1 1 A 1 

- = — + </>—■ 
r r„ 


(35) 


Another method has been to include photon recycling as another generation term in the 
continuity equation [79]. If we were to follow this approach, we could modify our 
continuity equation in Equation (17) to become 


^-u(r,t) = G dectron +G phalon -sd(vj) (36) 

OX beam recycling € 

Additional research would be necessary to determine the functional form of the photon 
recycling generation rate. 
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c. 


GENERATION VOLUME 


When an electron from an SEM enters a semiconductor, it undergoes a series of 
elastic and inelastic collisions that deflect it from its path and cause lateral spreading 
from its initial trajectory. The scattering events create a 3D generation volume of charge 
carriers below the surface on a pm scale [80]. Under high energy (e.g., 20 keV) electron 
beam excitation, the average energy required to make a single EHP is only 6.5 eV for 
TIBr and 4.64 eV for CZT. A single high energy electron can therefore create thousands 
of EHPs within the generation volume [21], [81]. At normal incidence of the electron 
beam, the shape of the generation volume depends on the incident beam energy and the 
atomic composition of the semiconductor. Monte Carlo electron simulations in low- 
density, low-Z materials show pear-shaped generation volumes, whereas simulations in 
high-density, high-Z materials reveal more hemispherical-shaped generation 
volumes [81], 

The 2D diffusion model for charge transport in bulk semiconductors models the 
generation volume as a delta function at a depth zo below the surface. When the diffusion 
lengths are large compared to the generation volume, the effects of the generation volume 
are minimal, and the delta function approximation is valid. When the diffusion lengths 
are comparable in size to the generation volume, however, the delta function 
approximation may fail to accurately model the generation volume. 

The net recorded luminescence in transport imaging is a convolution of the 
intensity distribution due to recombination in the generation volume and the intensity 
distribution due to charge transport. Luminescence from the generation volume depends 
on the energy of the incident electrons, whereas luminescence at points removed from the 
generation volume holds information about the diffusion behavior of the charge carriers. 
In order to de-convolve these two distributions, the luminescence distribution for the 
generation volume must be modeled. Two approaches considered here for future work 
are Monte Carlo simulations and experimentation. 

A Monte Carlo simulation code called CASINO, version 2.4.8.1, has been run to 
investigate the generation volume in bulk TIBr. The code models the trajectories and 
energy loss of primary electrons until they reach energies of 50 eV. Below 50 eV, the 
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modeling becomes inaccurate [82]. For the simulations in TIBr, the SEM beam energy is 
20 keV, the SEM spot size is 60 nm, and the number of simulated electrons is 100,000. 
Figure 65(a) shows the trajectory of primary electrons after scattering. Electrons in this 
simulation penetrate to a maximum depth of ~ 1700 nm. Figure 65(b) shows the 
distribution of the electrons as a function of their maximum depth. On average, the 
electrons reach a depth of ~700 nm, which was the value used for zo in all mathematical 
fitting for TIBr in this dissertation. Figure 65(c) shows a side view (XZ) of the energy 
density distribution in depth of the material integrated over the width of the material. 
Approximately 95%of the energy is deposited within ~ 750 nm of the surface. Figure 
65(d) shows a planar view (XY) of the lateral extent of the energy density distribution 
integrated through the depth of the material. Approximately ninety-five percent of the 
SEM beam energy is distributed in an area with a diameter of 
~500 pm. 



Figure 65. Monte Carlo simulations for 20 keV electrons in TIBr showing (a) Electron 
trajectories, (b) Electron depth distribution, (c) Integrated energy density 
distribution in XZ plane, and (d) Integrated energy density distribution in XY 

plane. 
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By plotting the energy density distributions shown in Figures 64(c) and (d) as 
functions of position, one can begin to create a model for the generation volume. Figure 
66 shows oblique and side views of the integrated energy density in TIBr. This energy 
density is the basis for the generation volume. 




Figure 66. Monte Carlo simulations for 20 keV electrons in TIBr showing (a) Oblique and 

(b) Side views of the integrated energy density. 


In Figure 67, the energy distribution shown in Figure 66 is fit to individual 
Gaussian and Lorentzian functions. Figure 67 indicates that a sum of Lorentzian or 
Gaussian functions may provide a good fit for the energy density. 


95 








Figure 67. Gaussian and Lorentzian fits to integrated energy distribution from Monte Carlo 

simulations with 20 keV electrons in TIBr. 

Monte Carlo simulations may provide a basis for future exploration into the 
generation volume. However, there are certain limitations to this approach. First, the 
Monte Carlo simulations cut off at 50 eV, thereby losing the full energy deposition for 
each electron. Second, the energy density distributions only reflect energy deposition, 
not the luminescence emitted from the recombination of electron-hole pairs created by 
this energy. The luminescence intensity distribution will be larger than the energy 
distribution. Much work is necessary to further model the luminescence emitted as a 
result of this energy deposition. Third, the energy density distributions are very narrow, 
with ~ 95% energy density deposition in an area with a diameter of - 100 nm, as shown 
in Figure 66(b). This distribution alone seems too narrow to account for our need to 
begin fitting the data at 4 pm due to the generation volume. 

An alternate experimental approach would be to use known transport imaging 
intensity distributions from high Z materials that just reflect the generation volume, with 
minimal contribution from carrier transport. Transport imaging results have been 
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reported on thick samples of highly doped / 2 -type GaAs at a probe current of 1><10' 9 A 
and beam energies of 5-30 keV, with the SEM beam operated in spot mode of excitation 
[81]. In that work, the small diffusion length of the material (< 1 pm) and the 3D nature 
of the carrier diffusion allowed for direct imaging of the generation volume. Figure 68(a) 
shows the results of that work. Figure 68(b) shows Gaussian and Lorentzian fits to the 
20 keV distribution. From Figure 68(b), the Lorentzian appears to be a reasonable fit to 
the generation volume. Interestingly, the generation volume distribution approaches zero 
at ~4 pm, the same distance used as a starting point for our mathematical fitting. This 
empirical distribution could be used for transport imaging in bulk samples to de-convolve 
the intensity distribution due to recombination in the generation volume and the intensity 
distribution due to diffusion of charge carriers in order to more accurately model the 
diffusion behavior of carriers. 



Figure 68. (a) Intensity distributions for n-type GaAs with beam energies 

from 5-30 kV, and (b) Gaussian and Lorentzian fits to 
the 20 keV intensity distribution. After [81]. 

D. HIGH-INJECTION MODEL 

Finally, the 2D diffusion model is based on the assumptions of low-injections 
conditions in a doped semiconductor, resulting in a constant carrier lifetime. Low 
injection conditions mean that the excess carrier concentration generated by the electron 
beam is much less than the equilibrium concentration, or n e « no. Doping allowed us to 
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say no = No- For TIBr at 5 K, however, the equilibrium carrier concentration is — 0 due to 
the large bandgap and low temperature, and n e » n 0 , resulting in high-injection 
conditions. We must therefore re-visit our assumptions of the model and consider how 
high injection affects carrier lifetime and recombination. One change to the model is 
considered here. 

The spontaneous radiative lifetime for band to band recombination was given 
previously as 


T — —— : 

R R 


NN 


Ak B T(n 0 +p 0 ) 


1 + 


«o + P> 


V 0 


0 J 


-1 


(30) 


Under low injection conditions where no = N D , Equation (30) reduced to 




BN, 


(34) 


Here, the lifetime is independent of the excess carrier concentration. If instead, however, 
we assume high injection conditions, where n e » no, and semi-insulating materials where 
no =po, Equation (30) reduces to 

1 


t r = 


Bn e( X ) ’ 


(35) 


Now the lifetime depends on the excess carrier concentration. Additionally, we have 
indicated that the excess carrier concentration has a spatial dependence (1-D case shown). 
The spatial dependence results from a rapid decrease in excess carrier concentration as 
the carrier concentration gradient drives diffusion away from the generation volume [37]. 
If we define the recombination rate as 

n„ (x) 


R = 


(36) 


then substitution of Equation (35) into Equation (36) results in a new expression for the 
recombination rate, given by 


R = B[n e (x)] 2 


(37) 


The recombination rate now has a quadratic dependence on the excess carrier 
concentration [38]. 
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The modified recombination rate is representative of one change to the model that 
should be considered for accurate modeling of semi-insulating materials in the high 
injection limit. This may be especially important in the case of a high-resistivity direct 
bandgap semiconductor, such as CZT. For TIBr, an indirect gap material where defect- 
mediated recombination is dominant, it appears likely that a constant lifetime model is 
appropriate, as indicated by the excellent fits obtained with the existing model for the Se- 
doped material. Further experimental studies as a function of doping would provide 
additional insight. 
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IX. CONCLUSION 


In this work, the transport properties of charge carriers in bulk semi-insulating 
semiconductors for nuclear radiation detectors have been investigated. Sample analysis 
has included undoped detector grade TIBr (45-S10) and doped TIBr (Na/Al/Ag, 
Cu/Fe/Zn, Se, and Pb), as well as CZT (ST1 and 372) under different post-growth 
conditions. Analytical techniques have included GDMS, EDX, TOF-SIMS, CL, PL, 
PICTS, and transport imaging using a 2D diffusion model for bulk materials. Using these 
methods we have demonstrated the following: 

• Low temperature CL spectrum in TIBr. 

• Application of a 2D model of carrier diffusion in bulk TIBr and CZT using 
transport imaging to advance a micro-characterization technique for rapid 
assessment of the magnitude and spatial variation in the jut product at 
2 pm resolution. 

• Spatial variation in the ambipolar diffusion length in TIBr on a scale of 
~10 pm. 

• Contribution to first empirical energy level diagram for defect levels of Se 
and Pb in TIBr. 

• Application of transport imaging in TIBr to demonstrate variations in 
charge transport properties in TIBr with doping/impurity levels. 

• Application of transport imaging in TIBr and CZT to report variations in 
/lit product with temperature over the ranges 8 K-102 K and 5 K-60 K, 
respectively. 

• Lirst application of transport imaging for preliminary investigation of the 
effects of Te inclusions on charge transport in CZT. 

• Lirst application of transport imaging for preliminary investigation of the 

presence of photon recycling in TIBr and CZT. 

101 



The magnitude and spatial variation of the /uz product are critical to the high 
resolution performance of room temperature semiconductor gamma-ray detectors. The 
research presented in this dissertation can contribute to the development and evaluation 
of these materials and devices by providing a rapid, contact-free, all-optical technique for 
micro-scale characterization of fundamental charge transport properties. 
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APPENDIX A. CL MAP 


The Matlab file CL_map.m opens a image file (.tif) and creates a CL map with a scale bar, as 
shown in Figure Al. 

CLmap.m 


clear 

clc 

% Load image file 
image = '0723 BSE'; 

10 = imread([image '.tif']); 

10 = double(10); % [X Width 2184, Y Height 1472] in MicroCCD 

sizelO = size (10); 

% Normalize 10 matrix 
m=max(max(10)); 

I0=I0/m; 

% Plot Map 
figure(1); 
surf(10); 
view([0,-90]); 
shading interp 
colormap(hot) 
colorbar 

xlim([0,sizeIO(2)]); ylim([0,sizeIO(l)]) 
set(gca, 'XTick' , [], 'XTickLabel' , [] ) 
set(gca, 'YTick' , [], 'YTickLabel' , [] ) 
saveas(gcf,[ 'CL_map_' image '-jpg']) 
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Figure Al. CL map. 
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APPENDIX B. CL SPECTRA 


The Matlab file CL_spectra.m opens one or more .dat files and creates one or more CL 
spectrum, as shown in Figures B2 and B3. The code shown here also includes Gaussian fits to 
several peaks. 

CL_spectra.m 


clc 

clear 

%% Load data file (.dat) and sample information 
fileload = { ... 

1, 'TIBr1' , 'TIBr 45-S10', '2000x', '6E-10A', '295K', '23FEB11'; 

2, 'TIBr9' , 'TIBr 45-S10', '2000x', '3E-10A', 'll K', '23FEB11'}; 

k = size(fileload,1); 
for j =1:k 

index{j} = fileload]j,1}; 

datafilefj} = [fileload{j,2} '.dat']; 

sample]j} = fileload]j,3}; 

mag{j} = fileload]j,4}; 

current{j} = fileload]j,5}; 

temp{j} = fileload]j,6}; 

date{j} = fileload]j,7}; 

end; 

%% Generate wavelength (X) and normalized intensity (Y) vectors 
for j = 1:k 

fid = fopen(datafile{j }); 

Header = textscan(fid, '%s' ,10, 'delimiter' ,' \n' ); 

Spectrum = textscan(fid, '%f %f'); 
fclose(fid); 

c = 0; % 0 = wavelength x axis, 1 = energy x axis 
if c == 0 

X{ j} = Spectrum]1,1}; 
xaxis = 'Wavelength (nm) '; 
elseif c == 1 

X{j} = 1.24 * 1000./Spectrum]1,1}; 
xaxis = 'Energy (keV) '; 

end 

Y{j} = Spectrum]1,2}; 

Ymax]j} = max(Y]j}); 

Ybkg]j} = min(Y{j}); 

Y{j} = (Y{j}-Ybkg]j})/(Yrnax]j}-Ybkg]j}); % Y values normalized to 
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end; 


%% Plot Spectrum 

xl=linspace(350,500,1001); 

xl0=416; 

al = 5.5; 

A1=0.84; 

fl=@(xl) Al*exp(-(abs(xl-xlO)/al). A 2); 

Gl=f1(xl); 

x2=linspace(350,500, 1001) ; 

x20=428; 

a2 = 10; 

A2 = 1; 

f2=@(x2) A2*exp(-(abs(x2-x20)/a2) . A 2) ; 

G2=f2(x2); 

x3=linspace(400,440,1001); 
x30=417; 
a3 = 10; 

A3=0.06; 

f3=0(x3) A3*exp(-(abs(x3-x30)/a3). A 2); 

G3=f3(x3); 

% TIBr(45-S10), CL 295K 
figure(1) 
hold on 

plot(X{l},Y{l}, 'b' , 'LineWidth' , 2 ) ; 
plot(xl,G1, '--r' , 'Markersize' ,2, 'LineWidth ',!); 
plot(x2 , G2 , '—g' , 'Markersize' ,2, 'LineWidth' , 1) ; 
hold off 

xlabel(xaxis, 'FontSize' ,20, 'FontWeight' , 'b' ) 

% xmin = 1.5; xmax = 3.6; xlim([xmin xmax]) 
xmin = 350; xmax = 800; xlim([xmin xmax]) 

ylabel( 'Intensity (Arb. units) ' , 'FontSize' , 20, 'FontWeight' , 'b' ) 

% ymin = -0.1; ymax = 1.1; ylim([ymin ymax]) 

set(gca, 'FontSize ',16, 'FontWeight ','b'); % set(gcf, 'Color', [1,1,1]) 

legend (' Data ',' Gaussian Fit 1', 'Gaussian Fit 2') 
set(gca, 'XMinorGrid' , 'off' , 'YGrid' , 'off' ) 
saveas(gcf, '295K.jpg') 

% TIBr (45-S10) , CL UK 
figure(2) 
hold on 

plot(X{2},Y{2}, 'b' , 'Markersize ' , 5, 'LineWidth' , 2) ; 
plot(x3,G3, '--r' , 'Markersize' , 2 , 'LineWidth' , 2) ; 
hold off 

xlabel(xaxis, 'FontSize' ,20, 'FontWeight' , 'b' ) 

% xmin = 1.5; xmax = 3.6; xlim([xmin xmax]) 
xmin = 300; xmax = 800; xlim([xmin xmax]) 

ylabel( 'Intensity (Arb. units )', 'FontSize' , 20, 'FontWeight' , 'b' ) 
ymin = -0.01; ymax = 1.1; ylim([ymin ymax]) 

set(gca, 'FontSize ',16, 'FontWeight ','b'); % set(gcf, 'Color', [1,1,1]) 
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Intensity (Arb. units) 


off' ) 


legend( 'Data' , 'Gaussian Fit 1 ) 

set(gca, 'XMinorGrid' , 'off', 'YGrid' , ' 

saveas(gcf, '11K.jpg' ) 



Figure Bl. CL spectrum at 11 K. 
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Figure B2. CL spectrum at 295 K. 
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APPENDIX C. INTENSITY DISTRIBUTION AND FIT 


This appendix explains how to create the intensity distribution plots and fit the model to the 
data for multiple files. 

Begin by loading and analyzing the image files. Use Loaddata_multifile.m. Enter the line image 
file names in the "File Information" section of Loaddata_multifile.m. 

The call structure is as follows: 

Loaddata_multifile.m 

Tiff2data4fit_multifile.m 

CenterCut_v2.m 

Norm4fit.m 

makeX.m 

sdvec.m 

reduce_outliars.m (optional) 

The result is a matlab file (.mat) with the filename "data_..." 

To plot the intensity distribution, use plot_distribution_RHS.m. This script loads and plots the 
data files, as shown in Figure Cl. 

To fit the data, use fit_LMA_LVA_vl_multifile.m. 

The call structure is as follows: 

fit_LMA_LVA_vl_multifile.m 

LMA_LVA_vl.m 

Ueval_v8.m 

UpartL_v7.m 

UpartV_v8.m 

polyEstimate.m 

Ueval_v8.m 

errorfit.m 

fitplot.m 

ConvergeCheck.m 

fitPlot.m 
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Loaddata multifile.m 


clc 
clear 
elf reset 

%% File Information 


fileload = { ... 

1, 'ST1 linel 180s 500x 6E-10A lOkV 8.12mm 295K x492.png', 

'CZT(STl)', 'Ll' ,~'x492 ' , '295K', T 22N0V11'; 

2, 'ST1 linel 180s 500x 6E-10A lOkV 8.12mm 295K x391.png', 

'CZT(STl)', 'Ll', 'x391', '295K', T 22N0V11'; 

3, 'ST1 linel 180s 500x 6E-10A lOkV 8.12mm 295K x732.png', 

'CZT(STl)', 'Ll' ,~'x732 ' , '295K', T 22N0V11'; %~OMIT~ 

4, 'ST1 line2 180s 500x 6E-10A lOkV 8.15mm 295K x322.png', 

'CZT(STl)', ~'L2' ,~'x322 ' , '295K', T 22N0V11'; 

5, 'ST1 line2 180s 500x 6E-10A lOkV 8.15mm 295K xl99.png', 

'CZT(STl)', 'L2', 'xl99', '295K', T 22N0V11'}; 

k = size(fileload,1); 


for j =1:k 

index{j}=fileload{j, 1}; 
imagefile{j}=fileload{j,2}; 
sample{j}=fileload{j,3}; 
location{j}=fileload{j,4}; 
line{j}=fileload{j,5}; 
temp{j}=fileload{j, 6}; 
date{j}=fileload{j,7}; 

datafile{j}=([ 'data ' sample{j} ' ' location{j} ' ' line{j} ' ' 

temp{j} ' ' date{j}]) ; 

end 


%% FUNCTION CALL FOR CUTTING OF THE IMAGE FILE AND SAVING DATA VECTOR 

for j = 1 % manually iterate through filenames for all files 

[DataVec X IO R varvec] = Tiff2data4fit multifile(imagefile{j}, 
sample{j}); 

save(datafile{j}, 'DataVec 'X','IO','R', 'varvec' ) 

end; 
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Tiff2data4fit multifile.m 


function [ DataVec X 10 R varvec] = Tiff2data4fit multifile( image, 
sample) 

10 = imread(image); 

10 = double(10); % [X Width 2184, Y Height 1472] in MicroCCD 

10 = 10'; 

%% Surface plot of data 
figurel = figure(1); 
elf( 'reset' ); 

axesl = axes(' Parent ',figurel ,' CLim ',[ 0 1]); 

view(axes 1, [0 0]) ; 

grid(axesl, 'minor' ); 

hold(axesl, 'all' ); 

surf1(10) 

shading interp 

colormap copper 

title(sample) 

xlabel( 'row' , 'VerticalAlignment' , 'cap' , 'HorizontalAlignment' , 'center' ); 
ylabel( 'column' , 'Visible' , 'off' , 'HorizontalAlignment' , 'center' ); 

%% First cut of Data 
clc 

% Rstartl = 0; 

% Rendl = 500; 

Rstartl = input (' \nlnput start range for desired Row data (1st cut): 

') ; 

Rendl = input (' Input end range for desired Row data: '); 
for I = 1:Rendl-Rstartl+1 

T1 (:,I) = 10(:,I+Rstartl); 

end 

10 = Tl; 

%% Surface plot of data 
elf( 'reset' ) ; 
figurel = figure(1); 

axesl = axes(' Parent ',figurel ,' CLim ',[ 0 1]); 

view(axesl, [0 0]); 

grid(axesl, 'minor' ); 

hold(axesl, 'all' ); 

surf1(10) 

shading interp 

colormap copper 

title(sample) 

xlabel( 'row' , 'VerticalAlignment' , 'cap' , 'HorizontalAlignment' , 'center' ); 
ylabel( 'column ', 'Visible', 'off', 'HorizontalAlignment' , 'center' ); 

%% Second cut of data in same direction 

Rstart2 = input (' \nlnput start range for desired Row data (2nd cut) : 
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Rend2 = input (' Input end range for desired Row data (2nd cut) : ') ; 

R=[Rstartl Rendl Rstart2 Rend2]; 

for I = 1:Rend2-Rstart2+1 

T2( : , I) = 10(:,I+Rstart2); 

end 

10 = T2; 

%% Surface plot of data 
elf( 'reset' ); 
figurel = figure(1); 

axesl = axes(' Parent figurel CLim ',[ 0 1]); 

view(axesl,[89 15]); 

grid(axesl, 'minor' ); 

hold(axesl, 'all' ); 

surf1(10) 

shading interp 

colormap copper 

title(sample) 

xlabel( 'row' , 'VerticalAlignment' , 'cap' , 'HorizontalAlignment' , 'center' ) 
ylabel( 1 column' , 'Visible', 'off', 'HorizontalAlignment ', 'center'); 

%% Save Ridge figure before CenterCut and X vec generation 

o, 

o 

% clf 

% figurel = figure (1); 

% axesl = axes('Parent',figurel,'CLim',[0 1]); 

% view(axesl, [0 0]) ; 

% grid(axesl,'minor'); 

% hold(axesl,'all'); 

% surf1(10) 

% shading interp 
% colormap copper 

% title(strcat(sample, ' (Before CenterCut)')) 

o, 

o 

xlabel('row','VerticalAlignment','cap','HorizontalAlignment','center') 
% ylabel('column','Visible','off','HorizontalAlignment','center'); 

O, 

o 

% %% Surface plot of data 
% % elf reset; 

% figure2 = figure (2); 

% axesl = axes('Parent',figure2,'CLim', [ 0 1]); 

% view(axesl, [89 15]); 

% grid(axesl,'minor'); 

% hold(axesl,'all'); 

% surf1(10) 

% shading interp 
% colormap copper 
% title(sample) 

o, 

o 

xlabel('row','VerticalAlignment','cap','HorizontalAlignment','center') 
% ylabel('column','Visible','off','HorizontalAlignment','center'); 
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%% Reducing outlayers 

% - if g is set to 0 then this section is skipped 

o, 

o 

% - if g is set to 1 then any value in that is more than 1 std away 
from 

% values in the row will be set to the mean of that row 

g = 0; 
if g == 1 

10 = reduce_outliars(10); 

end 


%% Centering and X vec generator 

Xmax = input( '\nEnter the number of pixels to the right of the peak to 
keep: ' ); 

% Xmax = 1000; 

[10] = CenterCut_v2(10,Xmax); 

% close(figure2); 
figure3 = figure (3); 

% elf('reset'); 

axesl = axes (' Parent figure3 CLim ',[ 0 1]); 

view(axesl,[90 0] ) ; 

grid(axesl, 'minor' ); 

hold(axesl, 'all' ); 

surf1(10) 

shading interp 

colormap copper 

title(sample) 

% title(strcat(sample, ' (After CenterCut)')) 

% title(strcat('Post-CenterCut Data for ', sample)) 

xlabel( 'row' , 'VerticalAlignment' , 'cap' , 'HorizontalAlignment' , 'center' ); 
ylabel( 1 column' , 'Visible' , 'off' , 'HorizontalAlignment' , 'center' ); 

%% Setting max value to 1 
10 = Norm4fit(10); 

%% Takes the mean of the rows of the matrix to create a single vector 
DataVec = mean(10'); 

%% Cuts the DataVec 

% so that the values of DataVec are above the noise and background 
levels 

[X] = makeX(DataVec); 
varvec = sdvec(X,10); 

%% Benediction 
display(' ' ) 

display( 'Goodbye!' ) 
display( ' ' ) 


113 



CenterCut v2.m 


function [ BB ] = CenterCut v2( AA , Xmax) 

% This will cut the matrix AA so that the columns will be of length 
Xmax 

% Position 1 of the AA matrix is the max value of that column 

[M,N]= size(AA); 

DV = mean(AA'); 
nl = max(DV) ; 

n2 = find((DV) == nl,1, 'first' ); 
for I = 1 : N 

n3 = find(AA(:,I) == max(AA(n2-20:n2+20,I)),1, 'last' ); 

BB(:,I) = AA(n3:n3+Xmax-l,I); 

end 

end 
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Norm4fit.m 


function [ BB ] = Norm4fit( AA ) 

% For each row in AA, a background noise level is calculated by taking 
a 

% mean of the furthest 20 value counts away from the source and then 
% subtracting this background noise level from each value in that row 

% The rows are then normalized to a max value of '1' 

[M, N] = size (AA) ; 

for I = M-19:M 

nl(I-M+20) = mean(AA(I,1:N)) ; 

end 

AA = AA-mean(nl); 

for 1=1:N 

n2=AA(1,I); 

BB(:,I)= AA(:,1) ./ (n2); 

end 

end 


makeX.m 

function [ X ] = makeX( DataVec ) 

% creates the X vector of same length as DataVec with .4 as stepsizes 

% N - length of DataVec 

[M,N] = size(DataVec); 

[M,N] = size(DataVec); 

X=0:N-1; 

X = X*(,4e-6); 

end 
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sdvec.m 


function [ varvec ] = sdvec( X ,10 ) 

[ml,nl] = size(X); 
varvec = zeros(ml,nl); 
for i=l:nl 

varvec(i)=std(10(i,:)); 

end 

end 


reduce_outliars.m (optional) 


function [ 10 ] = reduce_outliars( 10 ) 

%Any value in that is more than 2 std away from 
%values in the row will be set to the mean of that row 

for 1=1:size(10,1) 

a = 2*std(10(I,:)); 
b = mean(10(1,:)); 
for J=1:size (10,2) 
if 10(I,J)>=b+a 
10(I,J)=b; 

elseif I0(I,J)<=b-a 
10 (I, J) =b; 

end 

end 

end 

end 
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plot_distribution_RHS.m 


clc 
clear 
elf reset 

%% File Information 

fileload = { ... 

1, 'ST1 linel 180s 500x 6E- 

'CZT(STl)', 'Ll' ,~'x492 ' , '295K 

2, 'STl_linel_180s_500x_6E- 

'CZT(STl)', 'Ll', 'x391' , '295K 

3, 'STl_linel_180s_500x_6E- 

'CZT(STl)', 'Ll', 'x732' , '295K 

4, 'ST1 line2 180s 500x 6E- 

'CZT(STl)', ~'L2' ,~'x322 ' , '295K 

5, 'ST1 line2 180s 500x 6E- 

'CZT(STl)', 'L2', 'xl99', '295K 

k = size (fileload,1); 

for j =1:k 

index{j}=fileload{j, 1}; 
imagefile{j}=fileload{j,2}; 
sample{j}=fileload{j,3}; 
location{j}=fileload{j,4}; 
line{j}=fileload{j,5}; 
temp{j}=fileload{j, 6}; 
date{j}=fileload{j,7}; 
datafile{j}=([ 'data ' sample 
temp{j} ' ' date{j}]) ; 

end 

DATA = [index' sample' location' 
%% PLOT DISTRIBUTIONS 
for j = l:k 

S{j} = load(datafile{j}); 

X{ j} = S{j},X*le6; 

Y{j} = S{j}.DataVec; 

end; 

o, o, 
o o 

% 295K 

figure(1) 
semilogy( ... 

X{1},Y{1}, 'r' , . . . 
X{2},Y{2}, 'g' , . . . 

X{ 3 }, Y{ 3 }, 'b ' , . . . 
X{4},Y{4}, 'm' ,... 


0A lOkV 8.12mm 295K x492.png', 
, T 22N0V11' ; 

0A lOkV 8.12mm 295K x391.png', 
, t 22NOV11'; 

0A lOkV 8.12mm 295K x732.png', 
, t 22NOV11'; 

0A lOkV 8.15mm 295K x322.png', 
, t 22NOV11' ; 

0A lOkV 8.15mm 295K xl99.png', 
, t 22NOV11' }; 


j} ' ' location{j} ' ' line{j} ' ' 


line' temp' date'] 


1 

I 

1 

I 

1 

I 

1 

I 

1 

I 

{ 
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X { 5 }, Y { 5 }, ' k ' , . . . 

'Markersize' ,5, 'LineWidth' , 2 ) ; 
xlabel (' Position (\mum) ' , 'FontSize' , 20, 'FontWeight' , 'b' ) 
xmin = 0; xmax = 25; xlim([xmin xmax]) 

ylabel( 'Intensity (Arb. units)' , 'FontSize', 20, 'FontWeight' , 'b' ) 
ymin = le-3; ymax = 1; ylim([ymin ymax]) 

set(gca, 'FontSize ',16, 'fontWeight' , 'b' ); % set(gcf, 'Color', [1,1,1]); 

set(gca, 'XMinorGrid' , 'off' , 'YMinorGrid' , 'off' ) 


A = 

[sample]1} ' 

' location]1} ' 

' line{1} ' 

' temp{1}]; 

B = 

[sample{2] ' 

' location{2} ' 

' line{2} ' 

' temp{2}]; 

C = 

[sample]3} ' 

' location{3} ' 

' line{3} ' 

' temp{3}]; 

D = 

[sample]4} ' 

' location{4} ' 

' line{4} ' 

' temp{4}]; 

E = 

[sample]5} ' 

' location{5} ' 

' line{5} ' 

' temp{5}]; 


legend(A,B,C,D,E, 'Location', 'NorthEast' ); 
figurenamesave = [sample]1} ' ' sample{4} ' ' temp{l}]; 

saveas(gcf, [figurenamesave '.jpg']); 



Position (|im) 


Figure Cl. Intensity distribution. 
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fit LMA LVA vl multifile.m 


format short e 
format compact 

clc 

clear 

%% File Information 


fileload = { ... 

1, 'ST1 linel 180s 500x 6E-10A lOkV 8.12mm 295K x492.png', 

'CZT(STl)', 'Ll' ,~'x492 ' , '295K', T 22N0V11'; 

2, 'ST1 linel 180s 500x 6E-10A lOkV 8.12mm 295K x391.png', 

'CZT(STl)', 'Ll' , 'x391' , '295K', T 22N0V11'; 

3, 'ST1 linel 180s 500x 6E-10A lOkV 8.12mm 295K x732.png', 

'CZT(STl)', 'Ll', 'x732', '295K', T 22N0V11'; 

4, 'ST1 line2 180s 500x 6E-10A lOkV 8.15mm 295K x322.png', 

'CZT(STl)', ~'L2' ,~'x322 ' , '295K', T 22N0V11'; 

5, 'ST1 line2 180s 500x 6E-10A lOkV 8.15mm 295K xl99.png', 

'CZT(STl)', ~'L2' , 'xl99' , '295K', T 22N0V11'}; 


k = size(fileload,1); 


for j =1:k 

index{j}=fileload{j, 1}; 
imagefile{j}=fileload{j,2}; 
sample{j}=fileload{j,3}; 
location{j}=fileload{j,4}; 
line{j}=fileload{j,5}; 
temp{j}=fileload{j, 6}; 
date{j}=fileload{j,7}; 

datafile{j}=([ 'data ' sample{j} ' ' location{j} ' ' line{j} ' ' 

temp{j} ' ' date{j}]) ; 

end 


for j =1:k 

file_name_2_load{j} = datafile{j}; 

filename_4_save{j} = (['fit ' sample!j} ' ' location!j} ' ' line{j} 

' ' temp{j} ' ' date{j}]); 

filel{j} = load(file_name_2_load{j}); 

end 


%% INITIAL PARAMETERS 

% initial inputs for parameters are to be entered for pvec start vector 


L = 10; 

V = 8; 

Zo = 0.5; 
A = 40; 


% enter L in microns 
% enter V in 1/micron 
% enter Zo in microns 


o, o, 
o o 


##################################################################### 
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first data point 2 use = 11; % [11,15]=4-10 urn; [11,20] =4 — 12 

um; [11,30]=4-16 um; [ 11,40 ] =4-20 um; [11, 115] =4-50 urn; 
number of data points 2 use = 115; % [6, 15] =2—8 um; [6,25]=2-12; 

[6,35]=2-16; [6,45]=2-20 um; [6,50]=2-22 um; 


pvec = [L V Zo A] ; 
initial_pvec = pvec; 

%%##################################################################### 


for j =1:k 

Xvec[j} = filel{j}.X*le6; 
Dvec{j} = f ilel{j}.DataVec; 


cc=le-5; 

xvec = Xvec[j](first data point 2 use:first data point 2 use+ 
number of data points 2 use ) ; 

dvec = Dvec[j](first data point 2 use:first data point 2 use+ 
number of data points 2 use ) ; 

display( ' ' ); 

display(' Start - new loop'); 
display ( '- ' ) ; 

j j=i; 

[ L V A phi]=LMA LVA vl(pvec(1) , pvec(2) , pvec(3) , pvec(4) , xvec, 
dvec,cc ); 

pvec(1)=L; 
pvec(2)=V; 
pvec(3)=Zo; 
pvec(4)=A; 

[uest] = fitPlot(pvec(1), pvec(2), pvec(3), pvec(4), xvec, dvec, 
filename_4_save{j}); 

save(filename 4 save{j pvecfirst data point 2 use', 'number of data 
_points_2_use' , 'initial_pvec' , 1 xvec' , 'dvec' , 'uest' ) 

display( 'Goodbye' ) 


end 


%% Comments/notes 

% pvec - parameter vector [ L V Zo A] 

% M - matrix of most recent 10 parameters [phi ; L ; V ; Zo ; A] 
% phi - squarred error vector (evec.*evec) 

% evec - error vector 
% L - diffusion length in meters 
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o\° o\° o\° o\° o\° o\° o\° o\° o\° o\° o\° o\° o\° o\° o\° o\° o\° o\° o\° o\° o\° o\° o\° o\° o\° 


V - S/D in 1/sec 

Zo - penetration depth of electon beam for model 
A - Amplitude of model 


FIT - is the final matrix that contains the data for each of the fits 
Lfit - temporary vector to save corresponding final values for L 
Vfit - same - for V 

Afit - same for A 

phifit - same for phi 

Dvec - vector of with all data points recorded 
dvec - vector in use for fitting of current data set 

Xvec - position vector of all the x-values 

xvec - position vector in use for fitting current data set 


Uest - Full length vector of fitted values with current parameters 
uest - vector of fitted values cut to desird length 


UpL - vector of partial derrivatives with respect to L 
UpV - wrt V 
UpZ - wrt Z 
UpA - wrt A 


Ln 

Vn 

Zn 

An 

En 


normalized vector of partial derrivatives with respect to L 
wrt V 
wrt Z 
wrt A 

normalized error vector 


grad - gradient vector with the four paramters 
Ngrad - normalized gradient vector 


% first data_point 2 use - # of data points moving away from the source 
to skip 

% number of data_points 2 use - # of data points to be used in the 
fitting of the model 

% M - Matrix that contains the last 10 values throught the iteration 
% process for [phi pvec(l) pvec(2) pvec(3) pvec(4) phistep ] 

% cc - is the percentage of change in each paramter of pvec allowed to 
% determine convergence 

% gM - AMatrix of 4 colums consiting of the normalized partials 
derivatives 

% for UpL,UpV,UpZ,UpA 

% LMA - Levenberg-Marquardt-Algorithm step direction 


% stepvec - vector used to gage a the best next step between the max 
and 

% min of the step vector 
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o\° o\° o\° f+ o\° o\° h) o\° o\° o\° 


% stepVecVariance - the amount that will be used in varying phistep (as 
a 

percentage) when creating the next 'stepvec' 
stepvec start - used to keep track of the initial stepvec so that if 
max(stepvec) is greater than '1' stepvec will be 

set 

to 'stepvec_start' 

fitcount - the number of times that the fitting will be performed - 
e 

fit with the lowest phi value in all this will be used, 
fitmat - matrix used to store the fitted values for each iteration of 
fitcount 
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LMA LVA vl.m 


function [LVA phi b] = LMA LVA vl(L,V ,Zo,A ,xvec, dvec, cc) 
pvec=[L V Zo A]; 


format long e 

% pvec - parameter vector [ L V Zo A] 

% M - matrix of most recent 10 parameters [phi ; L ; V ; Zo ; A] 
% phi - squarred error vector (evec.*evec) 

% evec - error vector 


% Dvec - vector of with all data points recorded 
% dvec - vector in use for fitting of current data set 


% Xvec - position vector of all the x-values 

% xvec - position vector in use for fitting current data set 


% Uest - Full length vector of fitted values with current parameters 
% uest - vector of fitted values cut to desird length 


% UpL - vector of partial derrivatives with respect to 

% UpV - wrt V 

% UpZ - wrt Z 

% UpA - wrt A 


L 


% Ln 
% Vn 
% Zn 
% An 
% En 


normalized vector of partial derrivatives with respect to L 
wrt V 
wrt Z 
wrt A 

normalized error vector 


% grad - gradient vector with the four paramter 
% Ngrad - normalized gradient vector 

% dp2skip - # of data points moving away from the source to skip 
% dp2use - # of data points to be used in the fitting of the model 

% M - Matrix that contains the last 10 values throught the iteration 
% process for [phi pvec(l) pvec(2) pvec(3) pvec(4) phistep ] 

% cc - is the percentage of change in each paramter of pvec allowed to 
% determine convergence 

% gM - AMatrix of 4 colums consiting of the normalized partials 
derivatives 

% for UpL,UpV,UpZ,UpA 

% LMA - Levenberg-Marquardt-Algorithm step direction 

% stepvec - vector used to gage a the best next step between the max 
and 

% min of the step vector 

% stepVecVariance - the amount that will be used in varying phistep (as 
a 
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stepvec_start 


percentage) when creating the next 'stepvec' 

- used to keep track of the initial stepvec so that if 
max(stepvec) is greater than '1' stepvec will be 

reset 

% to 'stepvec_start' 

j j=0; 

converge=0; 

M = zeros(10,6) ; 
tic 

phistep = .9; 
stop=7; 
cc=le-6; 
see =5 ; 
lambda = leO; 
nu = 3 ; 
phi=0; 


while converge ==0; 

jj=jj+i; 

uest = Ueval_v8(pvec(1),pvec(2),pvec(3),pvec(4),xvec); 

UpL = UpartL_v7(pvec(1),pvec(2),pvec(3),pvec(4) , xvec); 

UpV = UpartV_v8(pvec(1),pvec(2),pvec(3),pvec(4), xvec); 

UpA = uest/pvec(4); 

evec = (uest-dvec); 

Aphi = evec*evec'; 

% al = [UpL ; UpV ; UpA ] ' ; 

al = [UpL/norm(UpL) ; UpV/norm(UpV); UpA/norm(UpA)]'; 
b = -(al'*al+lambda*diag(diag(al'*al)))\(al'*evec'); 

LMA = b'/norm(b); 

Bpvec = [pvec(1)+LMA(1)*phistep*pvec(1) 
pvec(2)+LMA(2)*phistep*pvec(2) pvec(3) pvec(4)+LMA(3)*phistep*pvec(4)] 
Buest = Ueval_v8(Bpvec(1),Bpvec(2),Bpvec(3),Bpvec(4),xvec); 

Bevec = Buest-dvec; 

Bphi = Bevec*Bevec'; 
phistepL=phistep; 
k=0 ; 

if BphikAphi 

if lambda <le-10 
Aphi=Bphi; 

end 

while Bphi<Aphi 
k=k+l; 

Aphi = Bphi; 

Blambda=lambda/(nu A k); 

b = -(al'*al+Blambda*diag(diag(al'*al)))\(al'*evec'); 

LMA = b'/norm(b); 

phistep = polyEstimate( pvec, LMA,xvec,dvec,phistepL,stop 
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Bpvec = [pvec(1)+LMA(1)*phistep*pvec(1) 
pvec(2)+LMA(2)*phistep*pvec(2) pvec(3) pvec(4)+LMA(3)*phistep*pvec(4)] 
Buest = Ueval_v8(Bpvec(1),Bpvec(2),Bpvec(3),Bpvec(4),xvec) 
Bevec = Buest-dvec; 

Bphi = Bevec*Bevec'; 
if Blambda <le-10 
Aphi=Bphi; 

end 

end 
k=k-l; 

else 

Clambda=lambda/(nu) ; 

b = -(al'*al+Clambda*diag(diag(al'*al)))\(al'*evec') ; 

LMA = b'/norm(b); 

phistep = polyEstimate( pvec, LMA,xvec,dvec,phistepL ,stop ); 
Cpvec = [pvec(1)+LMA(1)*phistep*pvec(1) 
pvec(2)+LMA(2)*phistep*pvec(2) pvec(3) pvec(4)+LMA(3)*phistep*pvec(4)] 
Cuest = Ueval_v8(Cpvec(1),Cpvec(2),Cpvec(3),Cpvec(4),xvec); 
Cevec = Cuest-dvec; 

Cphi = Cevec*Cevec'; 

if CphikAphi 
k=k+l; 

Bphi=Cphi; 
if lambda <le-10 
Aphi=Bphi; 

end 

while BphikAphi 
k=k+l; 

Aphi = Bphi; 

Blambda=lambda/(nu A k) ; 

b = -(al'*al+Blambda*diag(diag(al'*al)))\(al'*evec'); 
LMA = b'/norm(b); 

phistep = polyEstimate( pvec, LMA,xvec,dvec,phistepL 

,stop ); 

Bpvec = [pvec(1)+LMA(1)*phistep*pvec(1) 
pvec(2)+LMA(2)*phistep*pvec(2) pvec(3) pvec(4)+LMA(3)*phistep*pvec(4)] 
Buest = 

Ueval_v8 (Bpvec (1) , Bpvec (2) , Bpvec (3) , Bpvec (4) , xvec) ; 

Bevec = Buest-dvec; 

Bphi = Bevec*Bevec'; 
if Blambda <le-10 
Aphi=Bphi; 

end 

end 
k=k-l; 

else 

if lambda >lelO 
Aphi=Bphi; 

else 

k=k-l; 

end 

while AphikBphi 
k=k-l; 

Bphi=min(Aphi,Bphi); 
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Alambda=lambda/(nu A k); 

b = -(al'*al+Alambda*diag(diag(al'*al)))\(al'*evec'); 
LMA = b'/norm(b); 

phistep = polyEstimate( pvec, LMA,xvec,dvec,phistepL 

,stop ); 

Apvec = [pvec(1)+LMA(1)*phistep*pvec(1) 
pvec(2)+LMA(2)*phistep*pvec(2) pvec(3) pvec(4)+LMA(3)*phistep*pvec(4)] 
Auest = 

Ueval_v8(Apvec(1),Apvec(2),Apvec(3),Apvec(4) , xvec) ; 

Aevec = Auest-dvec; 

Aphi = Aevec*Aevec'; 
if Alambda >lelO 
Aphi=Bphi; 

end 

end 

Bpvec=Apvec; 
k=k+l; 

end 

end 

lambda=lambda/(nu A k) ; 

% al = [UpL/norm(UpL) ; UpV/norm(UpV); UpA/norm(UpA) ] 
b = -(al'*al + lambda*diag(diag(al'*al)))\(al'*evec') ; 

LMA = b'/norm(b); 

phistep = polyEstimate( pvec, LMA,xvec,dvec,phistepL ,stop ); 

pvec = [pvec(1)+LMA(1)*phistep*pvec(1) 
pvec(2)+LMA(2)*phistep*pvec(2) pvec(3) pvec(4)+LMA(3)*phistep*pvec(4)] 
uest = Ueval_v8(Bpvec(1),Bpvec(2),Bpvec(3),Bpvec (4),xvec); 
evec = uest-dvec; 
phi = evec*evec'; 


M (1:9,1:6)=M(2:10,1:6); 

M(10,:)=[phi pvec(l) pvec(2) pvec(3) pvec(4) phistep ]; 


converge = ConvergeCheck(M,cc); 

% if converge == 1 

% if stop <6 

% converge=0; 

% stop=stop+l; 

% cc=cc/5; 

% display(' ############################# YAHOOOO 

###########') 


o, 

o 

o, 

o 

Q. 

O 

O, 

O 

O, 

O 

O, 

o 

o, 

o 

o, 

o 

o, 

o 

o, 

o 


M = zeros(10,6); 

M(10,:)=[phi pvec(l) pvec(2) pvec(3) pvec(4) phistep ]; 

Lname = 'L = ' ; 

Vname = ', V = '; 

Zname = ', Zo = '; 

Aname = ', A = '; 

pi = 'step ='; 
p2 = ', phi ='; 
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% LI = ' , LMA = ' ; 

% lam = ', time ='; 

% laml = ', jj count ='; 

% lam2 = ', lambda ='; 

% timel= toe; 

% fprintf('%s %1.4e %s %1.5e %s %1.4e %s %1.4e %s %4.3f 

%10s %4.3e \n'... 

% ,Lname,pvec(1),Vname,pvec(2),Zname, pvec(3) , 

Aname,pvec(4), lam,timel,lam2,lambda') 

o, 

o 

% fprintf('%s %1.10f %s %1.9e %s %3.3e %s %4.4e %4.4e %4.4 

\n \n \n'... 


O, 

o 


o, 

o 


o, 

o 

O 

O, 

O 

o, 

o 

o, 

o 


LMA(3)); 
xvec,dvec. 


,pl, phistep,p2, phi,laml,jj, LI , LMA(l) 

fitPlot( pvec(l), pvec(2) ,pvec(3),pvec(4) 

test' ); 
pause(.01) 


end 

end 


LMA(2) 


if mod(jj,see)==1; 

Lname = 'L = ' ; 

Vname = ', V = ' ; 

Zname = ', Zo = ' ; 

Aname = ', A = ' ; 

pi = 'step =' ; 
p2 = ', phi =' ; 

LI = ', LMA = ' ; 
lam = ', time =' ; 
laml = ', jj count ='; 
lam2 = ', lambda ='; 
timel= toe; 

fprintf( '%s %1.4e %s %1.5e %s %1.4e %s %1.4e %s %4.3f %10s 
%4.3e \n' ... 

,Lname,pvec(1),Vname,pvec(2),Zname,pvec(3), Aname,pvec(4), 
lam,timel,lam2,lambda') 

fprintf( '%s %1.10f %s %1.9e %s %3.3e %s %4.4e %4.4e %4.4e \n 

\n \n' ... 

,pl, phistep,p2, phi,laml,jj, LI , LMA(l) , LMA(2) 

,LMA(3)); 

fitPlot( pvec(l), pvec(2) ,pvec(3),pvec(4) , xvec, dvec, 'test' ) 

pause(.01) 


end 


end 


L = pvec(1); 
V = pvec(2); 
Zo = pvec(3); 
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A = pvec(4); 


Lname = 'L = ' ; 

Vname = ' , V = ' ; 

Zname = ', Zo = ' ; 

Aname = ', A = ' ; 

pi = 'step =' ; 
p2 = ', phi =' ; 

LI = ', LMA = ' ; 
lam = ', time =' ; 
laml = ', jj count ='; 
lam2 = ', lambda ='; 
timel= toe; 

fprintf( '%s %1.4e %s %1.5e %s %1.4e %s %1.4e %s %4.3f %10s %4.3e 

\n' . . . 

,Lname,pvec(1),Vname,pvec(2),Zname,pvec(3), Aname,pvec(4), 
lam,timel,lam2,lambda') 

fprintf( '%s %1.10f %s %1.9e %s %3.3e %s %4.4e %4.4e %4.4e \n \n \n' . 

,pl, phistep,p2, phi,laml,jj, LI , LMA(1) , LMA(2) ,LMA(3)); 

fitPlot( pvec(l), pvec(2) ,pvec(3),pvec(4) ,xvec,dvec, 'test' ); 

pause (.01) 

end 
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Ueval v8.m 


function [ Uest ] = Ueval_v8( L,V,Zo,A,xvec ) 

% evaluates U by the provided inputs 

% this requires the evaulation of two integrals. 'f1' and 'f2' 
% ef (easy function) is used to siimplify fl and f2 

k = length(xvec); 

U1 = zeros(1,k); 

U2 = zeros(1,k); 

U3 = zeros(1,k); 

o, o, 
o o 


parfor n=l:k 
a=le-13; 
b=le-13; 

x = xvec(n); 

o, o, 
o o 

LL= x/L; 

ef = @(t) sqrt( t. A 2 - (x/L). A 2 ); 

fl = @(t) besselk(1,t).*ef(t); 

U1(n)= 

quadgk(f1, LL, Inf, 'RelTol',a, 'AbsTol' ,b, 'MaxIntervalCount' , 9999999); 


% f3 = @ (r) besselk (0, sqrt (r+(x/L) A 2) ).* (exp (V* (Zo-L. *sqrt (r) ))- 

1)./(2.*sqrt(r)); 

% LL= (Zo/L) A 2; 

% f3 = @(r) besselk(0,sqrt(r+x A 2)/L).*( exp(V*(Zo-sqrt(r)))-1 

) ./ (2*L.*sqrt(r)); 

% LL=Zo A 2; 

f3 = @(r) besselk(0,sqrt(r. A 2+x A 2)/L) .* ( exp(V*(Zo-r))-1 )./L; 
LL=Zo; 

U3(n)= 

quadgk(f3, LL, Inf, 'RelTol' , a, 'AbsTol' ,b, 'MaxIntervalCount' , 9999999); 
end 

Uest=(U1+U3)*A; 
end 
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UpartL_v7.m 


function [ UpL ] = UpartL v7( L,V,Zo,A,xvec ) 

% This function calculates the partial derrivative of 
% u(L,V,Z,A) with respect to L 

% takes as inputs the value for each variable L,V,Z,A 
% and a vector of concerened x values 

% sends as output a vector with the partial of u wrt L 
% evaluated at each of the x vaules in the xvec 
%display('UparL') 

%tic 

k = length(xvec); 

11 = zeros(1,k); 

12 = zeros(1,k); 

13 = zeros(1,k); 

a= 1 e -13 ; 
b=1e-13; 

parfor n=l:k 

x = xvec(n); 

LL= x/L; 

ef = @(t) sqrt( t. A 2 - (x/L) A 2 ); 

fl = @(t) besselk(2,t).*ef(t).*(x A 2/L A 3)./t; 

II(n)= 

quadgk(f1, LL, Inf, 'RelTol',a, 'AbsTol' ,b, 'MaxIntervalCount' , 9999999); 


% f2 = @ (r) besselk (0, sqrt (r+(x/L) A 2) ).* ( (- 

V* (r+(x/L) A 2) .*exp(V*(Zo-L.*sqrt(r)))./(r) ) + ... 

% (1- exp(V*(Zo-L.*sqrt(r))) ) .*x A 2./(L A 3*r. A (3/2)) )/2; 

% LL= (Z o/L) A 2; 


% f2 = @(r) (1-exp(V.*(Zo-sqrt(r))))./(2*L A 2.*sqrt(r)).*( 

besselk(1,sqrt(r+x A 2)./L).*(sqrt(r+x A 2)./L)-besselk(0,sqrt(r+x A 2)./L) 
) ; 

% LL= Z o A 2; 


f2 = @(r) 1/(L A 2)*(l-exp(V.*(Zo- 

r))).*(besselk(1,sqrt(r. A 2+x A 2)./(L)).*sqrt(r. A 2+x A 2)./(L)- 
besselk(0, sqrt(r. A 2+x A 2) ./(L) ) ) ; 

LL= Zo; 
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12(n)= 

quadgk(f2,LL,Inf, 'RelTol' ,a, 'AbsTol' ,b, 'MaxIntervalCount' ,9999999) 
end 

UpL=(11+12)*A; 
end 
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o\° o\° o\° o\° 


UpartV_v8.m 


function [ UpV ] = UpartV v8( L,V,Zo,A,xvec ) 

This function calculates the partial derrivative of 
u(L,V,Z,A) with respect to V 

takes as inputs the value for each variable L,V,Z,A 
and a vector of concerened x values 

% sends as output a vector with the partial of u wrt V 
% evaluated at each of the x vaules in the xvec 
%display('UpartV') 

%tic 

k = length(xvec); 

UpV = zeros(l,k); 


o, o, 
o o 


a=l00*eps; 
b=1e-2 3; 

parfor n=l:k 

intl=0; 
x = xvec(n); 

tt = 3; 

if tt == 1 

LL = (Zo/L) A 2; 

fl = @(r) besselk(0,sqrt(r+(x/L) A 2)).*(Zo- 
L.*sqrt(r)).*exp(V.*(Zo-L.*sqrt(r)))./(2.*sqrt(r)); 

elseif tt== 2 
LL = Z o A 2; 

fl = @(r) besselk(0,sqrt(r+x A 2)/L).*exp(V*(Zo-sqrt(r))).*(Zo- 
sqrt(r) ) ./(2*L.*sqrt(r)); 


else 


LL = Zo; 

fl = @(r) besselk(0,sqrt(r. A 2+x A 2)/L).*exp(V*(Zo-r)).*(Zo- 


r)./(L); 
end 


o. 


accl=10; 

ww=10; 

stop = min(LL*accl,ww); 
c = 500; 


o 


o 


o 


o. 


o 


o 


o. 


o 


o 


o 


i=l ; 

while LL<=stop; 


o 


inc = log(1+(i)/c) A 4; 
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% inti = inti + 

quadgk(fl,LL,LL+inc,'RelTol',a,'AbsTol',b,'MaxIntervalCount 9999999 ) 

O, 

o 

Q. 

O 

% i=i+l; 

% LL=LL+inc; 

% end 

UpV(n)= inti + 

quadgk(f1, LL, Inf, 'RelTol' ,a, 'AbsTol' ,b, 'MaxIntervalCount' , 9999999); 


end 

UpV=UpV*A; 

end 
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polyEstimate.m 


function [ phistep ] = polyEstimate( pvec, LMA,xvec,dvec ,a,bl) 

% uses the elements of'stepvec' to detemine the best stepsize in the 
% direction of 'LMA' for pvec. 

% returns the estimated stepsize 

% stepvec = [le-5 le-5+.05 le-5+.l ]; 

o, 

o 

% ercheck=zeros(1,length(stepvec)); 

o, 

o 

% for j=1:length(stepvec) 

% Pcheck=[pvec(1)+LMA(1)*stepvec(j)*pvec(1) 

pvec(2)+LMA(2)*stepvec(j)*pvec(2) pvec(3) 
pvec(4)+LMA(3)*stepvec(j)*pvec(4) ] ; 

% Ucheck=Ueval_v8(Pcheck(1),Pcheck(2),Pcheck(3),Pcheck(4),xvec); 

o, 

o 

% ercheck(j)=(dvec-Ucheck)*(dvec-Ucheck)'; 

% end 

% phistep=errorfit3(stepvec, ercheck) ; 

% end 

a=min(a*100, .08) ; 

%b=-7; 
b=-bl; 

stepvec = [10 A b 3*10 A b a/6 a/4 a/3 a]; 

%stepvec = [ 10 A b (10 A b)*1.845 max((a/219),(10 A b)*2.5743) 

max((a/111),(10 A b)*2.8) max((a/10),(10 A b)*3.9) max((a),(10 A b)*5) ] ; 

%stepvec = [10 A b 10 A b+a 10 A b+(a*2) 10 A b+(a*3) ]; 


ercheck=zeros(1,length(stepvec)); 

for j=l:length(stepvec) 

Pcheck=[pvec(1)+LMA(1)*stepvec(j)*pvec(1) 
pvec(2)+LMA(2)*stepvec(j)*pvec(2) pvec(3) 
pvec(4)+LMA(3)*stepvec(j)*pvec(4)]; 

Ucheck=Ueval v8(Pcheck(1),Pcheck(2),Pcheck(3),Pcheck(4),xvec); 

ercheck(j) = (dvec-Ucheck)*(dvec-Ucheck) ' ; 

end 

phistep=errorfit(stepvec, ercheck, b) ; 
end 
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errorfit.m 


function [t] = errorfit( a,f,b ) 
nl = length (f); 

S = zeros(nl,nl); 

S(:,1)=f; 

for i = 2:nl 

for j = 2:i 

S(i, j) = (S(i, j-l)-S(i-l,j-1))/(a(i)-a(i-j+l)) ; 

end 

end 

fl =@(x) S(l,l) +(x-a(1) ) .*S (2,2) + (x-a(1)) .* (x-a(2)) ,*S (3,3) +. . . 

(x-a(1)).*(x-a(2)).*(x-a(3)).*S(4,4)+(x-a(l)).*(x- 
a (2) ) . *(x-a (3)) . *(x-a (4)) .*S (5,5) ; 


x=min(a)/10 A b:max(a)/10 A b; 
x=x*10 A b; 

% plot(x,f1(x)) 

t=x(find(f1(x)==min(fl(x)),1, 'first' ) ) ; 


end 
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ConvergeCheck.m 


function [ converge ] = ConvergeCheck( M,cc ) 

% This function uses matrix M and convergence criteria 'c' to determine 
if 

% the new pvec paramters are with in tolerance for convergence 

% 'd' is used to set the previous itterations for checking note that a 
% value of 7 will check the rows 7-9 against the 10th row or in other 
% words the previous three iterations 

d=7 ; 

converge =0; 

if abs(max(M(d:9,2))-M(10,2))<=M(10,2)*cc 

if abs(max(M(d:9,3))-M(10,3))<=M(10,3)*cc 

if abs(max(M(d:9,4))-M(10,4))<=M(10,4)*cc 

if abs(min(M(d:9,2))-M(10,2))<=M(10,2)*cc 

if abs(min(M(d:9,3))-M(10,3))<=M(10,3)*cc 

if abs(min(M(d:9,4))-M(10,4))<=M(10,4)*cc 

converge = 1; 


end 

end 

end 


end 

end 

end 


end 
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fitPlot.m 


function [ uest ] = fitPlot( L, V ,Zo ,A ,xvec,dvec, filenamesave 
%UNTITLED2 Summary of this function goes here 
% Detailed explanation goes here 

acc=0; 

uest=Ueval_v8(L,V,Zo,A,xvec); 
evec = (uest-dvec); 
phi = evec*evec'; 

figure(20) 
subplot (2,1,1) 

plot(xvec,dvec, 'o', xvec,uest, 'MarkerSize' ,5) 
legend( 'data' , 'model' ) 

title(strcat(filenamesave, ' , 1st dp :', num2str(xvec(1)),' microns , 
dp: ' ,num2str(xvec(length(xvec))) , 'microns' ), 'FontSize' ,15); 
xlabel (' position (microns) ', 'FontSize ', 13) 
ylabel (' Intensity (a. u .)',' FontSize ', 13) 
d = max(uest); 
c= min(uest); 

text(xvec(floor(length(xvec)/2)),(c+(d-c)*9/10),[ ' L = ' 

, num2str(round(L*le2)*le-2), ' microns'] , 'FontSize ' , 15) ; 
text(xvec(floor(length(xvec) 12 )), (c+(d-c)*7/10),[ ' S/D = ' 

,num2str(round((V*100))/100), '/micron'] , 'FontSize' , 15) ; 
text(xvec(floor(length(xvec) 12 )) , (c+(d-c)*5/10),[ ' Zo = ' 

,num2str(round(Zo*le2)*le-2), ' microns' ], 'FontSize' ,15); 

text(xvec(floor(length(xvec)/2)),(c+(d-c)*3/10),[ ' Amplitude = ' 

,num2str(round(A*le4)*le-4) ] , 'FontSize' ,15); 

text(xvec(floor(length(xvec) / 6 )) ,(c+(d-c)*9/10),[ ' phi = ' 

,num2str(round(phi*le7)/le3), '1e — 4 ' ], ' FontSize' ,15) 

subplot (2,1,2) 

semilogy(xvec,dvec, 'o', xvec,uest, 'MarkerSize' ,5) 
legend( 'data' , 'model' ) 

xlabel(' position (microns) ', 'FontSize ', 13) 
ylabel (' Intensity (a. u. )',' FontSize ', 13) 
saveas(gcf,strcat(filenamesave, ' -jpg')) 
end 


last 
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Intensity (a. u.) Intensity (a. u.) 


fitCZT(STI) LI x492 295K 22NOV11 , 1st dp:4microns , last dp:50microns 



Figure C2. Mathematical fit. 
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APPENDIX D. SPATIAL DISTRIBUTION 


This appendix explains how to create an intensity distribution plot, divide the distribution into 
desired lengths (i.e., 5 pixel or 2 pm segments) and fit the model to the data for each segment. 

Begin by loading and analyzing each image file. Use Loaddata_multicut.m. Enter the line image 
file names in the "File Information" section of Loaddata_multicut.m. The number of pixels used 
is controlled by changing the value of the variable p in the file finecut.m. 

The call structure is as follows: 

Loaddata_multicut.m 

Tiff2data4fit_multicut.m 

CenterCut_v2.m 

Norm4fit.m 

finecut.m 

makeX.m 

The result is a matlab file (.mat) with the filename "dataMC..." 

To fit the data, use fit_LMA_LVA_vl_multicut.m. This file loads the "dataMC..." file, fits each 
distribution, and creates a file containing all the fitted diffusion lengths. This file is called "fit...L 
data.m". 

The call structure is as follows. 

fit_LMA_LVA_vl_multicut.m 

LMA_LVA_vl.m 

Ueval_v8.m 

UpartL_v7.m 

UpartV_v8.m 

polyEstimate.m 

Ueval_v8.m 

errorfit.m 

fitplot.m 

ConvergeCheck.m 

fitPlot.m 

To plot the diffusion length for each segment, use the file plot_Ld.m. Enter the file fit...L data.m 
in the "File Information" section. The file plot_Ld.m calls mu_tau.m. The resulting plots are 
shown in Figures D1 and D2. 
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Loaddata multicut.m 


clc 
clear 
elf reset 

%% File Information 


fileload = { ... 

1, 'ST1 linel 180s 500x 6E-10A lOkV 8.12mm 295K x492.png', 

'CZT(STl)', 'Ll ' , ' x4 92^' , '295K', T 22N0V11'; 

2, 'ST1 linel 180s 500x 6E-10A lOkV 8.12mm 295K x391.png', 

'CZT(STl)', 'Ll' , 'x391' , '295K', T 22N0V11'; 

3, 'ST1 linel 180s 500x 6E-10A lOkV 8.12mm 295K x732.png', 

'CZT(STl)', 'Ll ' ,~'x732 ' , '295K', T 22N0V11'; 

4, 'ST1 line2 180s 500x 6E-10A lOkV 8.15mm 295K x322.png', 

'CZT(STl)', ~'L2' ,~'x322 ' , '295K', T 22N0V11'; 

5, 'ST1 line2 180s 500x 6E-10A lOkV 8.15mm 295K xl99.png', 

'CZT(STl)', ~'L2' , 'xl99' , '295K', T 22N0V11'}; 

k = size (fileload,1); 


for j =1:k 

index{j}=fileload{j, 1}; 
imagefile{j}=fileload{j,2}; 
sample{j}=fileload{j,3}; 
location{j}=fileload{j,4}; 
line{j}=fileload{j,5}; 
temp{j}=fileload{j, 6}; 
date{j}=fileload{j,7}; 

datafile{j}=([ 'dataMC ' sample{j} ' ' location!j} ' ' line{j} ' ' 

temp{j} ' ' date{j}]) ; 

end 


%% Select file and analyze data 

j=i; 

[ DataVecMC X IO R ] = Tiff2data4fit multicut(imagefile{j}, sample!j}); 
save(datafile{j}, 'DataVecMC 'X', 'IO', 'R') 
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Tiff2data4fit multicut.m 


function [ DataVecMC X 10 R ] = Tiff2data4fit multicut( imagefile, 
sample ) 

10 = imread(imagefile); 

10 = double(10); % [X Width 2184, Y Height 1472] in MicroCCD 

10 = 10'; 

%% Surface plot of data 
figurel = figure(1); 
elf( 'reset' ); 

axesl = axes(' Parent figurel CLim ',[ 0 1]); 

view(axes 1, [0 0]) ; 

grid(axesl, 'minor' ); 

hold(axesl, 'all' ); 

surf1(10) 

shading interp 

colormap copper 

title(sample) 

xlabel( 'row' , 'VerticalAlignment' , 'cap' , 'HorizontalAlignment' , 'center' ); 
ylabel( 'column' , 'Visible' , 'off' , 'HorizontalAlignment' , 'center' ); 

%% First cut of Data 

% Rstartl = 0; 

% Rendl = 500; 

Rstartl = input (' \nlnput start range for desired Row data (1st cut): 

') ; 

Rendl = input (' Input end range for desired Row data (1st cut): '); 

for I = 1:Rendl-Rstartl+1 

T1 (:,I) = 10(:,I+Rstartl); 

end 

10 = Tl; 

% Surface plot of data 
elf(' reset' ) ; 
figurel = figure(1); 

axesl = axes (' Parent ', figurel ,' CLim ',[ 0 1]); 

view(axesl, [0 0]) ; 

grid(axesl, 'minor' ); 

hold(axesl, 'all' ); 

surf1(10) 

shading interp 

colormap copper 

title(sample) 

xlabel( 'row' , 'VerticalAlignment' , 'cap' , 'HorizontalAlignment' , 'center' ); 
ylabel( 'column ', 'Visible', 'off', 'HorizontalAlignment' , 'center' ); 

%% Second cut of data 

Rstart2 = input (' \nlnput start range for desired Row data (2nd cut): 
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Rend2 = input (' Input end range for desired Row data (2nd cut): ') ; 

R=[Rstartl Rendl Rstart2 Rend2]; 

for I = 1:Rend2-Rstart2+1 

T2( : , I) = 10(:,I+Rstart2); 

end 

10 = T2; 

% Surface plot of data 
elf( 'reset' ); 
figurel = figure(1); 

axesl = axes (' Parent figurel CLim ',[ 0 1]); 

view(axesl,[89 15]); 

grid(axesl, 'minor' ); 

hold(axesl, 'all' ); 

surf1(10) 

shading interp 

colormap copper 

title(sample) 

xlabel( 'row' , 'VerticalAlignment' , 'cap' , 'HorizontalAlignment' , 'center' ) 
ylabel( 'column' , 'Visible' , 'off' , 'HorizontalAlignment' , 'center' ); 

%% Centering 

Xmax = input( '\nEnter the number of pixels to the right of the peak to 
keep: ' ) ; 

% Xmax = 1000; 

[10] = CenterCut_v2(10,Xmax); 

elf(' reset' ); 
figurel = figure(1); 

axesl = axes( 'Parent figurel CLim ',[ 0 1 ] ) ; 

view(axesl, [0 0 ]); 

grid(axesl, 'minor' ); 

hold(axesl, 'all' ); 

surf1(10) 

shading interp 

colormap copper 

title(sample) 

xlabel( 'row' , 'VerticalAlignment' , 'cap' , 'HorizontalAlignment' , 'center' ) 
ylabel( 'column' , 'Visible' , 'off' , 'HorizontalAlignment' , 'center' ) ; 

pause 

%% Setting max value to 1 
10 = Norm4fit(10); 

elf( 'reset' ); 
figurel = figure(1); 

axesl = axes(' Parent ',figurel ,' CLim ',[ 0 1]); 
view(axesl, [0 0 ]); 
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grid(axesl, 'minor' ); 
hold(axesl, 'all' ); 
surf1(10) 
shading interp 
colormap copper 
title(sample) 

xlabel( 'row' , 'VerticalAlignment' , 'cap' , 'HorizontalAlignment' , 'center' ) 
ylabel( 'column' , 'Visible', 'off', 'HorizontalAlignment' , 'center' ) ; 

%% Average over 2 um 

[DataVecMC]=finecut(10); 

DataVecMC=DataVecMC'; 

%% Create X vector 

[X] = makeX(DataVecMC); 

%% Creat X vector 

% [X]=fineX(B); 

%% Benediction 
display( ' ' ) 

display( 'Goodbye!' ) 
display( ' ' ) 
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CenterCut v2.m 


function [ BB ] = CenterCut v2( AA , Xmax) 

% This will cut the matrix AA so that the columns will be of length 
Xmax 

% Position 1 of the AA matrix is the max value of that column 

[M,N]= size(AA); 

DV = mean(AA'); 
nl = max(DV) ; 

n2 = find((DV) == nl,1, 'first' ); 
for I = 1 : N 

n3 = find(AA(:,I) == max(AA(n2-20:n2+20,I)),1, 'last' ); 

BB(:,I) = AA(n3:n3+Xmax-l,I); 

end 

end 


144 




Norm4fit.m 


function [ BB ] = Norm4fit( AA ) 


% For each row in AA, a background noise level is calculated by taking 
a 

% mean of the furthest 20 value counts away from the source and then 
% subtracting this background noise level from each value in that row 

% The rows are then normalized to a max value of '1' 

[M, N] = size (AA) ; 

for I = M-19:M 

nl(I-M+20) = mean(AA(1,1:N)); 

end 

AA = AA-mean(nl); 

for 1=1:N 

n2=AA(1,I); 

BB(:,I)= AA(:,1)./(n2); 

end 

end 
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finecut.m 


function [ B ] = finecut( 10 ) 

% This function cuts averages the columns of the 10 matrix in sets of 
% 5 columns. Since each column represents 0.4 microns, this averages 
% the data over 2 microns. 

p=5; % pixel set size 

[m, n] = size (10) ; 

B=zeros(m,n/p); 

for i=l:m 

for k=l:n/p 

A(k)=mean(10(i,(k-1)*p+l:p*k)); 

B (i, k) = A (k) ; 

end 

end 

end 


makeX.m 

function [ X ] = makeX( DataVec ) 

% creates the X vector of same length as DataVec with .4 as stepsizes 

% N - length of DataVec 

[M,N] = size(DataVec); 

[M,N] = size(DataVec); 

X=0:N-1; 

X = X*(,4e-6); 

end 
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fit LMA LVA vl multicut.m 


format short e 
format compact 

clc 

clear 

%% File Information 


fileload = { ... 

1, 'ST1 linel 180s 500x 6E-10A lOkV 8.12mm 295K x492.png', 

'CZT(STl)', 'Ll' ,~'x492 ' , '295K', T 22N0V11'; 

2, 'ST1 linel 180s 500x 6E-10A lOkV 8.12mm 295K x391.png', 

'CZT(STl)', 'Ll', 'x391', '295K', T 22N0V11'; 

3, 'ST1 linel 180s 500x 6E-10A lOkV 8.12mm 295K x732.png', 

'CZT(STl)', 'Ll', 'x732', '295K', T 22N0V11'; 

4, 'ST1 line2 180s 500x 6E-10A lOkV 8.15mm 295K x322.png', 

'CZT(STl)', ~'L2' ,~'x322 ' , '295K', T 22N0V11'; 

5, 'ST1 line2 180s 500x 6E-10A lOkV 8.15mm 295K xl99.png', 

'CZT(STl)', ~'L2' , 'xl99' , '295K', T 22N0V11'}; 


k = size(fileload,1); 


for j =1:k 

index{j}=fileload{j, 1}; 
imagefile{j}=fileload{j,2}; 
sample{j}=fileload{j,3}; 
location{j}=fileload{j,4}; 
line{j}=fileload{j,5}; 
temp{j}=fileload{j, 6}; 
date{j}=fileload{j,7}; 

datafile{j}=([ 'dataMC ' sample{j} ' ' location!j} ' ' line{j} ' ' 

temp{j} ' ' date{j}]) ; 

end 


%% Select file and analyze data 

j=i; 

file name 2 load = datafile{j}; 
file = load(file name 2 load); 

%% INITIAL PARAMETERS 

% initial inputs for parameters are to be entered for pvec start vector 

L = 10; % enter L in microns 

V = 8; % enter V in 1/micron 

Zo = 0.5; % enter Zo in microns 

A = 40; 

first data point 2 use = 11; 

number of data points 2 use = 115; % (11-40=4-20 urn) (11-115=4-40 urn) 
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pvec = [L V Zo A] ; 
initial_pvec = pvec; 


n=20; % Must match the number of rows in the DataVecMC vector 


for k=l:n 


Xvec = file.X*le6; 

DataVec(k,:) = file.DataVecMC(k,:); 

Dvec = DataVec(k,:); 

cc=le-5; 

xvec = Xvec(first data point 2 use:first data point 2 use+ 
number of data points 2 use); 

dvec = Dvec(first data point 2 use:first data point 2 use+ 
number of data points 2 use); 

display( ' ' ); 

display(' Start - new loop'); 
display ( '- ' ) ; 

j j=i; 

[ L V A phi]=LMA LVA vl(pvec(1),pvec(2),pvec(3),pvec(4), xvec, 
dvec,cc ); 

pvec(1)=L; 
pvec(2)=V; 
pvec(3)=Zo; 
pvec(4)=A; 

Lvec(k)=L; 

Yvec(k)=2*k-l; 

Diffvec(k,1) = Yvec(k); 

Diffvec(k,2) = Lvec(k); 

filename_4_save = (['fit ' sample{j} ' ' location{j} ' ' temp{j} ' 

' date{j} ' cut ' num2str(k)]); 

fitPlot( pvec(l), pvec(2), pvec(3), pvec(4), xvec, dvec, 
filename 4 save ); 


save(filename 4 save,'pvec','first data point 2 use','number of data po 
ints 2 use' , 'initial_pvec' ) 

display( 'Goodbye' ) 


end 
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save (['fit L data ' samplefj} ' ' location{j} ' 

date{j}], 'Diffvec' ) 


' temp{j} ' ' 
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LMA LVA vl.m 


function [LVA phi b] = LMA LVA vl(L,V ,Zo,A ,xvec, dvec, cc) 
pvec=[L V Zo A]; 


format long e 

% pvec - parameter vector [ L V Zo A] 

% M - matrix of most recent 10 parameters [phi ; L ; V ; Zo ; A] 
% phi - squarred error vector (evec.*evec) 

% evec - error vector 


% Dvec - vector of with all data points recorded 
% dvec - vector in use for fitting of current data set 


% Xvec - position vector of all the x-values 

% xvec - position vector in use for fitting current data set 


% Uest - Full length vector of fitted values with current parameters 
% uest - vector of fitted values cut to desird length 


% UpL - vector of partial derrivatives with respect to 

% UpV - wrt V 

% UpZ - wrt Z 

% UpA - wrt A 


L 


% Ln 
% Vn 
% Zn 
% An 
% En 


normalized vector of partial derrivatives with respect to L 
wrt V 
wrt Z 
wrt A 

normalized error vector 


% grad - gradient vector with the four paramter 
% Ngrad - normalized gradient vector 

% dp2skip - # of data points moving away from the source to skip 
% dp2use - # of data points to be used in the fitting of the model 

% M - Matrix that contains the last 10 values throught the iteration 
% process for [phi pvec(l) pvec(2) pvec(3) pvec(4) phistep ] 

% cc - is the percentage of change in each paramter of pvec allowed to 
% determine convergence 

% gM - AMatrix of 4 colums consiting of the normalized partials 
derivatives 

% for UpL,UpV,UpZ,UpA 

% LMA - Levenberg-Marquardt-Algorithm step direction 

% stepvec - vector used to gage a the best next step between the max 
and 

% min of the step vector 

% stepVecVariance - the amount that will be used in varying phistep (as 
a 
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stepvec_start 


percentage) when creating the next 'stepvec' 

- used to keep track of the initial stepvec so that if 
max(stepvec) is greater than '1' stepvec will be 

reset 

% to 'stepvec_start' 

j j=0; 

converge=0; 

M = zeros(10,6) ; 
tic 

phistep = .9; 
stop=7; 
cc=le-6; 
see =5 ; 
lambda = leO; 
nu = 3 ; 
phi=0; 


while converge ==0; 

jj=jj+i; 

uest = Ueval_v8(pvec(1),pvec(2),pvec(3),pvec(4),xvec); 

UpL = UpartL_v7(pvec(1),pvec(2),pvec(3),pvec(4) , xvec); 

UpV = UpartV_v8(pvec(1),pvec(2),pvec(3),pvec(4), xvec); 

UpA = uest/pvec(4); 

evec = (uest-dvec); 

Aphi = evec*evec'; 

% al = [UpL ; UpV ; UpA ] ' ; 

al = [UpL/norm(UpL) ; UpV/norm(UpV); UpA/norm(UpA)]'; 
b = -(al'*al+lambda*diag(diag(al'*al)))\(al'*evec'); 

LMA = b'/norm(b); 

Bpvec = [pvec(1)+LMA(1)*phistep*pvec(1) 
pvec(2)+LMA(2)*phistep*pvec(2) pvec(3) pvec(4)+LMA(3)*phistep*pvec(4)] 
Buest = Ueval_v8(Bpvec(1),Bpvec(2),Bpvec(3),Bpvec(4),xvec); 

Bevec = Buest-dvec; 

Bphi = Bevec*Bevec'; 
phistepL=phistep; 
k=0 ; 

if BphikAphi 

if lambda <le-10 
Aphi=Bphi; 

end 

while Bphi<Aphi 
k=k+l; 

Aphi = Bphi; 

Blambda=lambda/(nu A k); 

b = -(al'*al+Blambda*diag(diag(al'*al)))\(al'*evec'); 

LMA = b'/norm(b); 

phistep = polyEstimate( pvec, LMA,xvec,dvec,phistepL,stop 
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Bpvec = [pvec(1)+LMA(1)*phistep*pvec(1) 
pvec(2)+LMA(2)*phistep*pvec(2) pvec(3) pvec(4)+LMA(3)*phistep*pvec(4)] 
Buest = Ueval_v8(Bpvec(1),Bpvec(2),Bpvec(3),Bpvec(4),xvec) 
Bevec = Buest-dvec; 

Bphi = Bevec*Bevec'; 
if Blambda <le-10 
Aphi=Bphi; 

end 

end 
k=k-l; 

else 

Clambda=lambda/(nu) ; 

b = -(al'*al+Clambda*diag(diag(al'*al)))\(al'*evec') ; 

LMA = b'/norm(b); 

phistep = polyEstimate( pvec, LMA,xvec,dvec,phistepL ,stop ); 
Cpvec = [pvec(1)+LMA(1)*phistep*pvec(1) 
pvec(2)+LMA(2)*phistep*pvec(2) pvec(3) pvec(4)+LMA(3)*phistep*pvec(4)] 
Cuest = Ueval_v8(Cpvec(1),Cpvec(2),Cpvec(3),Cpvec(4),xvec); 
Cevec = Cuest-dvec; 

Cphi = Cevec*Cevec'; 

if CphikAphi 
k=k+l; 

Bphi=Cphi; 
if lambda <le-10 
Aphi=Bphi; 

end 

while BphikAphi 
k=k+l; 

Aphi = Bphi; 

Blambda=lambda/(nu A k) ; 

b = -(al'*al+Blambda*diag(diag(al'*al)))\(al'*evec'); 
LMA = b'/norm(b); 

phistep = polyEstimate( pvec, LMA,xvec,dvec,phistepL 

,stop ); 

Bpvec = [pvec(1)+LMA(1)*phistep*pvec(1) 
pvec(2)+LMA(2)*phistep*pvec(2) pvec(3) pvec(4)+LMA(3)*phistep*pvec(4)] 
Buest = 

Ueval_v8 (Bpvec (1) , Bpvec (2) , Bpvec (3) , Bpvec (4) , xvec) ; 

Bevec = Buest-dvec; 

Bphi = Bevec*Bevec'; 
if Blambda <le-10 
Aphi=Bphi; 

end 

end 
k=k-l; 

else 

if lambda >lelO 
Aphi=Bphi; 

else 

k=k-l; 

end 

while AphikBphi 
k=k-l; 

Bphi=min(Aphi,Bphi); 
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Alambda=lambda/(nu A k); 

b = -(al'*al+Alambda*diag(diag(al'*al)))\(al'*evec'); 
LMA = b'/norm(b); 

phistep = polyEstimate( pvec, LMA,xvec,dvec,phistepL 

,stop ); 

Apvec = [pvec(1)+LMA(1)*phistep*pvec(1) 
pvec(2)+LMA(2)*phistep*pvec(2) pvec(3) pvec(4)+LMA(3)*phistep*pvec(4)] 
Auest = 

Ueval_v8(Apvec(1),Apvec(2),Apvec(3),Apvec(4) , xvec) ; 

Aevec = Auest-dvec; 

Aphi = Aevec*Aevec'; 
if Alambda >lelO 
Aphi=Bphi; 

end 

end 

Bpvec=Apvec; 
k=k+l; 

end 

end 

lambda=lambda/(nu A k) ; 

% al = [UpL/norm(UpL) ; UpV/norm(UpV); UpA/norm(UpA) ] 
b = -(al'*al + lambda*diag(diag(al'*al)))\(al'*evec') ; 

LMA = b'/norm(b); 

phistep = polyEstimate( pvec, LMA,xvec,dvec,phistepL ,stop ); 

pvec = [pvec(1)+LMA(1)*phistep*pvec(1) 
pvec(2)+LMA(2)*phistep*pvec(2) pvec(3) pvec(4)+LMA(3)*phistep*pvec(4)] 
uest = Ueval_v8(Bpvec(1),Bpvec(2),Bpvec(3),Bpvec (4),xvec); 
evec = uest-dvec; 
phi = evec*evec'; 


M (1:9,1:6)=M(2:10,1:6); 

M(10,:)=[phi pvec(l) pvec(2) pvec(3) pvec(4) phistep ]; 


converge = ConvergeCheck(M,cc); 

% if converge == 1 

% if stop <6 

% converge=0; 

% stop=stop+l; 

% cc=cc/5; 

% display(' ############################# YAHOOOO 

###########') 


o, 

o 

o, 

o 

Q. 

O 

O, 

O 

O, 

O 

O, 

o 

o, 

o 

o, 

o 

o, 

o 

o, 

o 


M = zeros(10,6); 

M(10,:)=[phi pvec(l) pvec(2) pvec(3) pvec(4) phistep ]; 

Lname = 'L = ' ; 

Vname = ', V = '; 

Zname = ', Zo = '; 

Aname = ', A = '; 

pi = 'step ='; 
p2 = ', phi ='; 
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% LI = ' , LMA = ' ; 

% lam = ', time ='; 

% laml = ', jj count ='; 

% lam2 = ', lambda ='; 

% timel= toe; 

% fprintf('%s %1.4e %s %1.5e %s %1.4e %s %1.4e %s %4.3f 

%10s %4.3e \n'... 

% ,Lname,pvec(1),Vname,pvec(2),Zname, pvec(3) , 

Aname,pvec(4), lam,timel,lam2,lambda') 

o, 

o 

% fprintf('%s %1.10f %s %1.9e %s %3.3e %s %4.4e %4.4e %4.4 

\n \n \n'... 


O, 

o 


o, 

o 


o, 

o 

O 

O, 

O 

o, 

o 

o, 

o 


LMA(3)); 
xvec,dvec. 


,pl, phistep,p2, phi,laml,jj, LI , LMA(l) 

fitPlot( pvec(l), pvec(2) ,pvec(3),pvec(4) 

test' ); 
pause(.01) 


end 

end 


LMA(2) 


if mod(jj,see)==1; 

Lname = 'L = ' ; 

Vname = ', V = ' ; 

Zname = ', Zo = ' ; 

Aname = ', A = ' ; 

pi = 'step =' ; 
p2 = ', phi =' ; 

LI = ', LMA = ' ; 
lam = ', time =' ; 
laml = ', jj count ='; 
lam2 = ', lambda ='; 
timel= toe; 

fprintf( '%s %1.4e %s %1.5e %s %1.4e %s %1.4e %s %4.3f %10s 
%4.3e \n' ... 

,Lname,pvec(1),Vname,pvec(2),Zname,pvec(3), Aname,pvec(4), 
lam,timel,lam2,lambda') 

fprintf( '%s %1.10f %s %1.9e %s %3.3e %s %4.4e %4.4e %4.4e \n 

\n \n' ... 

,pl, phistep,p2, phi,laml,jj, LI , LMA(l) , LMA(2) 

,LMA(3)); 

fitPlot( pvec(l), pvec(2) ,pvec(3),pvec(4) , xvec, dvec, 'test' ) 

pause(.01) 


end 


end 


L = pvec(1); 
V = pvec(2); 
Zo = pvec(3); 
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A = pvec(4); 


Lname = 'L = ' ; 

Vname = ' , V = ' ; 

Zname = ', Zo = ' ; 

Aname = ', A = ' ; 

pi = 'step =' ; 
p2 = ', phi =' ; 

LI = ', LMA = ' ; 
lam = ', time =' ; 
laml = ', jj count ='; 
lam2 = ', lambda ='; 
timel= toe; 

fprintf( '%s %1.4e %s %1.5e %s %1.4e %s %1.4e %s %4.3f %10s %4.3e 

\n' . . . 

,Lname,pvec(1),Vname,pvec(2),Zname,pvec(3), Aname,pvec(4), 
lam,timel,lam2,lambda') 

fprintf( '%s %1.10f %s %1.9e %s %3.3e %s %4.4e %4.4e %4.4e \n \n \n' . 

,pl, phistep,p2, phi,laml,jj, LI , LMA(1) , LMA(2) ,LMA(3)); 

fitPlot( pvec(l), pvec(2) ,pvec(3),pvec(4) ,xvec,dvec, 'test' ); 

pause (.01) 

end 
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Ueval v8.m 


function [ Uest ] = Ueval_v8( L,V,Zo,A,xvec ) 

% evaluates U by the provided inputs 

% this requires the evaulation of two integrals. 'f1' and 'f2' 
% ef (easy function) is used to siimplify fl and f2 

k = length(xvec); 

U1 = zeros(1,k); 

U2 = zeros(1,k); 

U3 = zeros(1,k); 

o, o, 
o o 


parfor n=l:k 
a=le-13; 
b=le-13; 

x = xvec(n); 

o, o, 
o o 

LL= x/L; 

ef = @(t) sqrt( t. A 2 - (x/L). A 2 ); 

fl = @(t) besselk(1,t).*ef(t); 

U1(n)= 

quadgk(f1, LL, Inf, 'RelTol',a, 'AbsTol' ,b, 'MaxIntervalCount' , 9999999); 


% f3 = @ (r) besselk (0, sqrt (r+(x/L) A 2) ).* (exp (V* (Zo-L. *sqrt (r) ))- 

1)./(2.*sqrt(r)); 

% LL= (Zo/L) A 2; 

% f3 = @(r) besselk(0,sqrt(r+x A 2)/L).*( exp(V*(Zo-sqrt(r)))-1 

) ./ (2*L.*sqrt(r)); 

% LL=Zo A 2; 

f3 = @(r) besselk(0,sqrt(r. A 2+x A 2)/L) .* ( exp(V*(Zo-r))-1 )./L; 
LL=Zo; 

U3(n)= 

quadgk(f3, LL, Inf, 'RelTol' , a, 'AbsTol' ,b, 'MaxIntervalCount' , 9999999 ); 
end 

Uest=(U1+U3)*A; 
end 
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UpartL_v7.m 


function [ UpL ] = UpartL v7( L,V,Zo,A,xvec ) 

% This function calculates the partial derrivative of 
% u(L,V,Z,A) with respect to L 

% takes as inputs the value for each variable L,V,Z,A 
% and a vector of concerened x values 

% sends as output a vector with the partial of u wrt L 
% evaluated at each of the x vaules in the xvec 
%display('UparL') 

%tic 

k = length(xvec); 

11 = zeros(1,k); 

12 = zeros(1,k); 

13 = zeros(1,k); 

a= 1 e -13 ; 
b=1e-13; 

parfor n=l:k 

x = xvec(n); 

LL= x/L; 

ef = @(t) sqrt( t. A 2 - (x/L) A 2 ); 

fl = @(t) besselk(2,t).*ef(t).*(x A 2/L A 3)./t; 

II(n)= 

quadgk(f1, LL, Inf, 'RelTol',a, 'AbsTol' ,b, 'MaxIntervalCount' , 9999999); 


% f2 = @ (r) besselk (0, sqrt (r+(x/L) A 2) ).* ( (- 

V* (r+(x/L) A 2) .*exp(V*(Zo-L.*sqrt(r)))./(r) ) + ... 

% (1- exp(V*(Zo-L.*sqrt(r))) ) .*x A 2./(L A 3*r. A (3/2)) )/2; 

% LL= (Z o/L) A 2; 


% f2 = @(r) (1-exp(V.*(Zo-sqrt(r))))./(2*L A 2.*sqrt(r)).*( 

besselk(1,sqrt(r+x A 2)./L).*(sqrt(r+x A 2)./L)-besselk(0,sqrt(r+x A 2)./L) 
) ; 

% LL= Z o A 2; 


f2 = @(r) 1/(L A 2)*(l-exp(V.*(Zo- 

r))).*(besselk(1,sqrt(r. A 2+x A 2)./(L)).*sqrt(r. A 2+x A 2)./(L)- 
besselk(0, sqrt(r. A 2+x A 2) ./(L) ) ) ; 

LL= Zo; 
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12(n)= 

quadgk(f2,LL,Inf, 'RelTol' ,a, 'AbsTol' ,b, 'MaxIntervalCount' ,9999999) 
end 

UpL=(11+12)*A; 
end 
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o\° o\° o\° o\° 


UpartV_v8.m 


function [ UpV ] = UpartV v8( L,V,Zo,A,xvec ) 

This function calculates the partial derrivative of 
u(L,V,Z,A) with respect to V 

takes as inputs the value for each variable L,V,Z,A 
and a vector of concerened x values 

% sends as output a vector with the partial of u wrt V 
% evaluated at each of the x vaules in the xvec 
%display('UpartV') 

%tic 

k = length(xvec); 

UpV = zeros(l,k); 


o, o, 
o o 


a=l00*eps; 
b=1e-2 3; 

parfor n=l:k 

intl=0; 
x = xvec(n); 

tt = 3; 

if tt == 1 

LL = (Zo/L) A 2; 

fl = @(r) besselk(0,sqrt(r+(x/L) A 2)).*(Zo- 
L.*sqrt(r)).*exp(V.*(Zo-L.*sqrt(r)))./(2.*sqrt(r)); 

elseif tt== 2 
LL = Z o A 2; 

fl = @(r) besselk(0,sqrt(r+x A 2)/L).*exp(V*(Zo-sqrt(r))).*(Zo- 
sqrt(r) ) ./(2*L.*sqrt(r)); 


else 


LL = Zo; 

fl = @(r) besselk(0,sqrt(r. A 2+x A 2)/L).*exp(V*(Zo-r)).*(Zo- 


r)./(L); 
end 


o. 


accl=10; 

ww=10; 

stop = min(LL*accl,ww); 
c = 500; 


o 


o 


o 


o. 


o 


o 


o. 


o 


o 


o 


i=l ; 

while LL<=stop; 


o 


inc = log(1+(i)/c) A 4; 
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% inti = inti + 

quadgk(fl,LL,LL+inc,'RelTol',a,'AbsTol',b,'MaxIntervalCount 9999999 ) 

O, 

o 

Q. 

O 

% i=i+l; 

% LL=LL+inc; 

% end 

UpV(n)= inti + 

quadgk(f1, LL, Inf, 'RelTol' ,a, 'AbsTol' ,b, 'MaxIntervalCount' , 9999999); 


end 

UpV=UpV*A; 

end 
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polyEstimate.m 


function [ phistep ] = polyEstimate( pvec, LMA,xvec,dvec ,a,bl) 

% uses the elements of'stepvec' to detemine the best stepsize in the 
% direction of 'LMA' for pvec. 

% returns the estimated stepsize 

% stepvec = [le-5 le-5+.05 le-5+.l ]; 

o, 

o 

% ercheck=zeros(1,length(stepvec)); 

o, 

o 

% for j=1:length(stepvec) 

% Pcheck=[pvec(1)+LMA(1)*stepvec(j)*pvec(1) 

pvec(2)+LMA(2)*stepvec(j)*pvec(2) pvec(3) 
pvec(4)+LMA(3)*stepvec(j)*pvec(4) ] ; 

% Ucheck=Ueval_v8(Pcheck(1),Pcheck(2),Pcheck(3),Pcheck(4),xvec); 

o, 

o 

% ercheck(j)=(dvec-Ucheck)*(dvec-Ucheck)'; 

% end 

% phistep=errorfit3(stepvec, ercheck) ; 

% end 

a=min(a*100, .08) ; 

%b=-7; 
b=-bl; 

stepvec = [10 A b 3*10 A b a/6 a/4 a/3 a]; 

%stepvec = [ 10 A b (10 A b)*1.845 max((a/219),(10 A b)*2.5743) 

max((a/111),(10 A b)*2.8) max((a/10),(10 A b)*3.9) max((a),(10 A b)*5) ] ; 

%stepvec = [10 A b 10 A b+a 10 A b+(a*2) 10 A b+(a*3) ]; 


ercheck=zeros(1,length(stepvec)); 

for j=l:length(stepvec) 

Pcheck=[pvec(1)+LMA(1)*stepvec(j)*pvec(1) 
pvec(2)+LMA(2)*stepvec(j)*pvec(2) pvec(3) 
pvec(4)+LMA(3)*stepvec(j)*pvec(4)]; 

Ucheck=Ueval v8(Pcheck(1),Pcheck(2),Pcheck(3),Pcheck(4),xvec); 

ercheck(j) = (dvec-Ucheck)*(dvec-Ucheck) ' ; 

end 

phistep=errorfit(stepvec, ercheck, b) ; 
end 
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errorfit.m 


function [t] = errorfit( a,f,b ) 
nl = length (f); 

S = zeros(nl,nl); 

S(:,1)=f; 

for i = 2:nl 

for j = 2:i 

S(i, j) = (S(i, j-l)-S(i-l,j-1))/(a(i)-a(i-j+l)) ; 

end 

end 

fl =@(x) S(l,l) +(x-a(1) ) .*S (2,2) + (x-a(1)) .* (x-a(2)) ,*S (3,3) +. . . 

(x-a(1)).*(x-a(2)).*(x-a(3)).*S(4,4)+(x-a(l)).*(x- 
a (2) ) . *(x-a (3)) . *(x-a (4)) .*S (5,5) ; 


x=min(a)/10 A b:max(a)/10 A b; 
x=x*10 A b; 

% plot(x,f1(x)) 

t=x(find(f1(x)==min(fl(x)),1, 'first' ) ) ; 


end 
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ConvergeCheck.m 


function [ converge ] = ConvergeCheck( M,cc ) 

% This function uses matrix M and convergence criteria 'c' to determine 
if 

% the new pvec paramters are with in tolerance for convergence 

% 'd' is used to set the previous itterations for checking note that a 
% value of 7 will check the rows 7-9 against the 10th row or in other 
% words the previous three iterations 

d=7 ; 

converge =0; 

if abs(max(M(d:9,2))-M(10,2))<=M(10,2)*cc 

if abs(max(M(d:9,3))-M(10,3))<=M(10,3)*cc 

if abs(max(M(d:9,4))-M(10,4))<=M(10,4)*cc 

if abs(min(M(d:9,2))-M(10,2))<=M(10,2)*cc 

if abs(min(M(d:9,3))-M(10,3))<=M(10,3)*cc 

if abs(min(M(d:9,4))-M(10,4))<=M(10,4)*cc 

converge = 1; 


end 

end 

end 


end 

end 

end 


end 
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fitPlot.m 


function [ uest ] = fitPlot( L, V ,Zo ,A ,xvec,dvec, filenamesave 
%UNTITLED2 Summary of this function goes here 
% Detailed explanation goes here 

acc=0; 

uest=Ueval_v8(L,V,Zo,A,xvec); 
evec = (uest-dvec); 
phi = evec*evec'; 

figure(20) 
subplot (2,1,1) 

plot(xvec,dvec, 'o', xvec,uest, 'MarkerSize' ,5) 
legend( 'data' , 'model' ) 

title(strcat(filenamesave, ' , 1st dp :', num2str(xvec(1)),' microns , 
dp: ' ,num2str(xvec(length(xvec))) , 'microns' ), 'FontSize' ,15); 
xlabel (' position (microns) ', 'FontSize ', 13) 
ylabel (' Intensity (a. u .)',' FontSize ', 13) 
d = max(uest); 
c= min(uest); 

text(xvec(floor(length(xvec)/2)),(c+(d-c)*9/10),[ ' L = ' 

, num2str(round(L*le2)*le-2), ' microns'] , 'FontSize ' , 15) ; 
text(xvec(floor(length(xvec) 12 )), (c+(d-c)*7/10),[ ' S/D = ' 

,num2str(round((V*100))/100), '/micron'] , 'FontSize' , 15) ; 
text(xvec(floor(length(xvec) 12 )) , (c+(d-c)*5/10),[ ' Zo = ' 

,num2str(round(Zo*le2)*le-2), ' microns' ], 'FontSize' ,15); 

text(xvec(floor(length(xvec)/2)),(c+(d-c)*3/10),[ ' Amplitude = ' 

,num2str(round(A*le4)*le-4) ] , 'FontSize' ,15); 

text(xvec(floor(length(xvec) / 6 )) ,(c+(d-c)*9/10),[ ' phi = ' 

,num2str(round(phi*le7)/le3), '1e — 4 ' ], ' FontSize' ,15) 

subplot (2,1,2) 

semilogy(xvec,dvec, 'o', xvec,uest, 'MarkerSize' ,5) 
legend( 'data' , 'model' ) 

xlabel(' position (microns) ', 'FontSize ', 13) 
ylabel (' Intensity (a. u. )',' FontSize ', 13) 
saveas(gcf,strcat(filenamesave, ' -jpg')) 
end 


last 
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plot_Ld.m 


clc 

clear 

% format bank 

%% Load Diffusion Length Data 

load( 'fit L data CZT(STl) LI 295K 22NOV11') 
Y=Diffvec( : , 1) ; 

L=Diffvec( : , 2 ) ; 

k=size(L,1) ; 

for j = 1: k 

MT(j) = mu_tau(L (j),5) ; 

end; 

MT=MT'; 

%% Plot Diffusion Length vs. Position 


% Diffusion Length 

figure(l); % Common Marker Types: o,s,d,p,h 

plot(Y, L, ' :ob' , 'Markersize ',10, 'MarkerEdgeColor' , 'b' , 'MarkerFaceColor' , 

' b' ) ; 

xmin = 0; xmax = 40; xlim([xmin xmax]) 
ymin = 0; ymax = 10; ylim([ymin ymax]) 

xlabel (' Position (\mum)' , 'FontSize ',20, 'FontWeight' , 'b' ) 

ylabel (' Diffusion Length (\mum)FontSize ', 20, 'FontWeightb' ) 

set(gca, 'XTick' ,0:2:40) 

set(gca, 'Xgrid' , 'on'); 

set (gca, 'Ygrid' , 'on'); 

set(gca, 'FontSize ',16, 'FontWeight' , 'b ' ) ; 
legend( 'CZT(ST1) ' ) 
fitsave = (['plot L' ]) ; 
saveas(gcf, [fitsave '.jpg']); 

% Mobility-Lifetime Product 

figure(2); % Common Marker Types: o,+,*,.,s,d,p,h 

plot(Y,MT, ' :or' , 'Markersize' , 10, 'MarkerEdgeColor' , 'r' , 'MarkerFaceColor' 

, ' r ' ) ; 

xlabel (' Position (\mum)' , 'FontSize ',20, 'FontWeight' , 'b' ) 

ylabel (' \mu\tau Product (cm A 2/ V) ',' FontSize ', 20, 'FontWeight ','b ' ) 

set(gca, 'XTick' ,0:2:40) 

set (gca, 'Xgrid' , 'on' ); 

set (gca, 'Ygrid' , 'on'); 

set(gca, 'FontSize ',16, 'FontWeight' , 'b' ); 
legend( 'CZT(ST1)' ) 
fitsave = (['plot MT'] ); 
saveas(gcf, [fitsave '.jpg']); 

% figure(3) 
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plotyy(Y,L,Y,MT) 


% figure(4) 

% [AX,HI,H2 ] = plotyy (Y, L, Y,MT, 'plot'); 

% set(get (AX( 1),'Ylabel'),'StringDiffusion Length 
(\mum) ', 'FontSize' , 20, 'FontWeight' , 'b') 

% set(get (AX (2 ),' Ylabel'),'String\mu\tau Product (cm A 2/ 
V) ', 'FontSize',20, 'FontWeight' , 'b') 


xlabel('Position (\mum)', 

'FontSize' 

,20, 

'FontWeight', 

'b 

set(AX(1), 

'XTick',0:2:40, 

'FontSize' 

kO 
\ —1 

'FontWeight', 

'b 

set(AX(2), 

'XTick',0:2:40, 

'FontSize' 

, 16, 

'FontWeight', 

'b 


% % set(AX(1),'YTick',0:0. 1 :3, 'FontSize', 16, 'FontWeightb') 

% % set(AX(2), 'YTick',0:0.0004:0.0048, 'FontSize',16, 'FontWeight', 'b') 

% % set(AX(1),'YLim',[1.6 2.4]) 

% % set(AX(2),'YLim',[5.2E-5 11.2E-5]) 

O, 

o 

set(HI, 'LineStyleMarker', 'o', 'MarkerSize',8, 'MarkerEdgeColor', 'b 
', 'MarkerFaceColor' , 'b') 

o, 

o 

set(H2, 'LineStyleMarker' , 's', 'MarkerSize' , 8, 'MarkerEdgeColor', 'g 
', 'MarkerFaceColor' , 'g') 

% set(gca,'Xgrid','on'); 

% fitsave = ( [ 'plotyy_example' ] ); 

% saveas(gcf, [fitsave '.jpg']); 


mu tau.m 


function [ MT ] = mu tau(L,T) 

% The function mu-tau(L,T) calculates the mobility-lifetime product. 

% It accepts the diffusion length L in urn and temperature T in K as its 
% arguments and passes back the mobility-lifetime product MT 

k boltzman = 8.617e-5; % eV/K 

MT = (L A 2/(k boltzman*T))*le-8; 

end 
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Figure Dl. Plot of diffusion length as a function of position. 
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Figure D2. Plot of /u t product as a function of position. 
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APPENDIX E. MONTE CARLO SIMULATIONS 


Monte Carlo modeling was performed using CASINO v.2.4.8.1, with the following 
parameters: 

Sample Definition : TIBr, 1x10 nm thickness, no substrate, density 7.56 g/cm 

Microscope Setup : 20 keV beam, 1x10 s electrons simulated, 50 nm beam radium 

Distributions : Kanaya Okayama max range parameters, X Divisions: 250, Y Divisions: 

250, Z Divisions: 250 

Options : Default used 

Change Physical Models : Defaults used 

Once the simulation is complete, click on the “Export” button and save the data as a .dat 
file. The file may be opened with Notebook and saved in Excel. 

The data may be analyzed in Matlab using the fde read_CASINO.m. The fde creates a 
3x3x3 matrix, with the order of (z,x,y). There is a test tile, called Test2_xyz.txt, that 
may be run to understand how Matlab stores and displays the CASINO data. There is 
also an annotated version of Test2_xyz in Excel called Test2_xyz.xlsx that also shows 
how Matlab stores and displays the CASINO data. 

To run the simulation, enter the filename in read_CASINO.m under the section “Load 
CASINO Energy Data.” 

read CASINO.m 


clc 

clear 

%% Load CASINO Energy Data 

% Read XZ planes and build 3x3 CASINO matrix 

% filename = 'Test2 xyz.txt'; 

% x = 5; y = 4; z = 6; 

filename = 'TIBr 20keV lE5nm lE5H.dat'; 
x = 250; y = 250; z = 250; 

CASINO = ones(z,x,y); 

a = 11; % Row containing XZ Plane 0, where cell A1 = (0,0) 

for n = 1:y 

RS = a+(2*n)+z*(n-1); RE = RS+(z-l); 

XZ = dlmread(filename, '\t ',[RS,1,RE,x]); 

CASINO(:,:,n) = XZ; 
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disp([num2str(n) '_done']) 

end 

% Sum through z-direction 
AA = ones(y,x); 
for j = l:y 

XZ = CASINO ( :, : , j ) ; 

AA(j,:) = sum(XZ,l); 

end 

% Find indices of maximum value in matrix 
[maxA,ind] = max(AA(:)); 

[R,C] = ind2sub(size(AA),ind) ; 

% Normalize AA matrix 
M = max(max(AA)); 

AA = AA/M; 

save([filename '.mat']) 

% %% Load CASINO CL data 

o, 

o 

% filename2 = 'CL by_R.dat'; 

o, 

o 

% b = 2; % Row containing Radius = 0, where cell A1 = (0,0) 

O, 

o 

% Rstart = b; Rend = Rstartf (x-1); 

o, 

o 

% CL = dlmread(filename2, '\t[Rstart,0,Rend,1]); 

o, 

o 

% radius LR = CL(:,1); % Creates column vector of position (radius) 
data 

% radius RL = -flipud(radius LR); % Flips position vector 
% radius = [radius_LR; radius_RL]; % Concatenates radius vectors 

o, 

o 

% offset = x/2; 

% radius = radius+offset; 

o, 

o 

% intensity LR = CL(:,2); % Creates column vector of CL intensity data 
% intensity^LR = intensity LR/(max(intensity LR)); % Normalizes CL 
intensity vector 

% intensity RL = flipud(intensity^LR); % Flips position vector 
% intensity = [intensity LR; intensity RL]; % Concatenates radius 
vectors 

%% Plot data 

figure (1) 

X = -x/2:1:x/2-1; 

Y = -y/2:1:y/2-l; 
surf(X,Y,AA) 

% view([0 90]) 

% view( [ 0 0]) 
view([-37.5,30]) 


% XY 
% XZ 

% oblique 
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shading interp 

xlabel( 'Position (nm) ' , 'FontSize ',20, 'FontWeight' , 'b' ) 

ylabel( 'Position (nm) ' , 'FontSize ',20, 'FontWeight' , 'b' ) 

zlabel (' Energy Density (Arb. units )',' FontSize ', 20 ,' FontWeightb' ) 

set(gca, ' FontSize ',16, 'FontWeight' , 'b' ) 

% colorbar 

xlim([-x/2 x/2]); ylim([-x/2 x/2]) 

% saveas(gca,'Planar.jpg') 
saveas(gca, 'Oblique.jpg' ) 


figure (2) 

BB = AA (R, : ) ; 

X = -x/2:1:x/2-1; 

% X = 1:1: x ; 

plot(X,BB, ' — .b' , 'MarkerSize' , 4 ) 
legend (' Energy Density') 
xlim([min(X) max(X)]) 

% figure(3) 

% plot(radius, intensity,'sk','MarkerSize',4) 

% legend('CL Intensity (%)') 

% xlim([0 300]) 

% figure(4) 

% hold on 

% plot(X,BB,'ob','MarkerSize',4); legend('Energy Density') 

% plot(radius, intensity, 'sk','MarkerSize',4); legend('CL Intensity 

(%) ' ) 

% xlim([0 300]) 
figure (5) 

mesh(CASINO(:,:,R)); 
view([0,-90]) 

% xlim([35 65]); ylim([0 25]) 
figure(6) 

contourf(CASINO(:,:,R),6); 
view([0,-90]) 

xlim([100 200]); ylim([0 100]) 

% xlim([35 65]); ylim([0 25]) 

% figure(5) 

% [C,h] = contour(CASINO(:,:,R),6); 

% clabel(C,h) 

% view([0,-90]) 

% xlim([100 200]); ylim([0 100]) 


syms r 
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r = -x/2:1:x/2-1; 

% r = 0:0.1:x; 
rO = 0; 

% Lorentzian 
g = 3; 

L = g./((r-rO) . A 2+g A 2) ; 

Lmax = max(L); 

L = L/Lmax; 

% figure(5) 

% plot(r,L,' r ' ) ; 

% legend('Lorentzian'); 

% Gaussian 
a = 2 ; 

G = exp(-((abs((r-rO)/a A 2)) , A 2)) ; 

% figure(6) 

% plot(r, G, 'g' ) ; 

% title('Gaussian'); 

% Comparison 
figure(7) 
hold on 

plot(X,BB, 'ob' , 'MarkerSize' ,5, 'LineWidth' ,2) 
plot(r,G, 'g', 'LineWidth' , 2) ; 
plot(r,L, 'r', 'LineWidth' , 2 ) ; 

% plot(radius,intensity,'sk','MarkerSize',4) 

xlabel( 'Position (nm)' , 'FontSize ',20, 'FontWeight' , 'b' ) 

ylabel (' Energy Density (Arb. units) ') 

set(gca,' FontSize' ,16, 'FontWeight' , 'b' ) 

legend( 'Energy Density' , 'Gaussian' , 'Lorentzian' ) 

% legend('Energy Density', 'Gaussian', 'Lorentzian', 'CL Intensity (%) ') 
xlim([-125 125]) 
hold off 

saveas(gca, 'Fit.jpg' ) 


% plane = 100; 

Q, 

O 

% figure(1) 

% surf(CASINO(:,:,plane)); 

% view([180,90]) 

% figure (2) 

% contour(CASINO(:,:,plane)); 

% view([180,90]) 

o, 

o 

% figure(3) 

% [C,h] = contour(CASINO(:,:,plane),10); 
% clabel(C,h) 

% view([180,90]) 

o, 

o 

% figure (4); 
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% offset = max(CASINO(:)); 

% hold on; 

% for i = 1:50 : y 

% surf(CASINOi) + offset*{i—1)) 

% end 


Two sample output graphs are shown below. 



Position (nm) 


-100 


-100 


Position (nm) 


Figure El. Oblique view. 



Figure E2. Side view with Gaussian and Lorentzian fits. 
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